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16552 SETTLING VELOCITY OF SPHERICAL PARTICLES INCALM 
KEY WORDS: ‘Drag; Experimental data; ; Particles; Sediment transport; 
ABSTRACT: The settling velocities of round and smooth spherical particles in calm — 
water at 200 /LT R /LT 6,000 were analysed utilizing modern measuring techniques. 
When the present results are compared with the “standard” CTS = f(R) curve, ol 
- fall on the average 20 percent higher. However, the “standard” curve was established 
with the majority of data derived from fixed spheres tested in wind tunnels. Because > 
the experiments are limited, a new drag coefficient cannot be made. To settle a 
- existing discrepancy between “fixed” and “falling” spheres, it will be necessary to do 
_ further research, Possibly focusing on the effect of the Strouhal number on the drag of 


REFERENCE: Boillat, Jean (Research Asst., Lab. Hydraulique, Ecole 
Polytechnique Federal, CH-1015 Lausanne, Switzerland), and Graf, Walter H., 
‘Settling Velocity of Spherical Particles in Calm Water,” Journal of the Hydraulics 
_ Division, ASCE, Vol. 107, No. HY10, Proc. Paper | 16552, October, 1981, pp. 1123- 


H LENGTH OF FLOW SEPARATION OVER DUNES _ 


| KEY WORDS: Bed roughness; Dimensional analysis; Dunes; Dune sands; __ 
' _ Flow separation; Flow visualization; Grain size; Open channel flow; Sands; 
| Two dimensional flow; Uniform flow _ 
ABSTRACT: Dimensional analysis was used to define the ada of the flow saparation 
for the case of two-dimensional, uniform flow over dunes, which was expressed _ 
1 terms of bed-form geometry and the sand size of the bed material. Experiments were 
_ conducted in a rectangular flume using artificial dunes and a flow visualization method | 
to measure the flow separation length. When dunes are flat, length of the flow 
_ separation is affected by the sand-grain roughness. When dunes are steep, the length of =—s_—© 
_ the flow separation is due mainly to the bed-form geometry and is relatively constant. _ 
_ The reverse rectangular steps can be considered as dunes of infinite length, and a b 
- length of the flow separation for s such steps will always be poy than that obtained 


| REFERENCE: Engel, Peter (Research ‘Bags, Hydr. Sect., Hydr. Div., 
| National Water Research Inst., Canada Center for Inland Waters., P.O. Box 5050, 5 
Burlington, Ontario, Canada L7R 4A6), “Length of Flow Separation over Dunes,” 


{| Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY10, Proc. Paper 16549, 


16541 ALGORITHMS FOR PIPE NETWORK ANALYSIS 
WORDS: Algorithms; Convergence; Converging» flow; Network - 
analysis; Network analyzers; Piping systems; Reliability; Water distribution = 
ABSTRACT: Algorithms for analyzing steady-state flow conditions in pipe networks 
are described for general applications. The algorithms are based on both loop equations 
expressed in terms of unknown flowrates and node equations expressed in terms of : 
_ unknown heads. Five methods, which represent those in significant use today, are 
_ included. Using an extensive data base, describing a variety of pipe networks, the 
reliabilities of these algorithms for pipe network analysis were investigated. Numerous © 
_ convergence and reliability problems were documented. Two methods based on loop 
_ equations have superior convergence characteristics. Methods based on node equations 
less reliable; single adjustment methods must be employed with caution. 
REFERENCE: Wood, Don Prof... Dept. of Civ Engrg... Univ. of Kentucky, 
_ Lexington, Ky. 40506), and Rayes, A. G., “Reliability of Algorithms for Pipe Network 
_ Analysis,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY10, | Proc. 


7 
4 
| 
| 


| 
KEY WORDS: Boundary conditions; Calculations; Contraction; Discharge 


(water); Flow characteristics; Laplace equation; Numerical analysis; Sluice 
ABSTRACT: The classical problems of steady-state flow under a sluice gate and over a 
spillway are solved by the boundary integral equation method. In both cases an 
iterative methods is used to adjust the free-surface profiles until they satisfy both the 
kinematic and dynamic boundary conditions. In the case of sluice gate flow, a point by 
point adjustment method is used which quickly converges to the correct solution 
within three to six iterations. The results closely check previous, less efficient | 
calculations. Some of the variation is past experimental results is explained by | 
calculations in which the geometry of the gate is altered slightly near the tip. Even a 
‘small change from the ideal geometry can have a major effect on the contraction © 
coefficient. Flow over a spillway is a more difficult problem due to poor convergence — 
properties in the free-surface profile behind te 


‘REFERENCE: Cheng, Alexander H-D. (Grad. Student, School of Civ. and 

_ Environmental Engrg., Cornell Univ., Ithaca, N.Y. 14853), Liggett, James A., and Liu, _ 
Philip L-F., “Boundary Calculations of Sluice and Spillway Flows,” Journal of the — 
Hydraulics Division, Vol. 107, No. HY10, Proc. Peper 16540, October, 1981, 
pp. 1163-1178 


16535 JET INJECTIONS FOR OPTIMUM MIXING IN PIPE FLOW oF 
| WORDS: Diffusion; Discharge measurement; Discharge (water); 
— Jet flow; Mixing; Optimization; Pipe flow; Secondary flow; Solutes 
ABSTRACT: The use of jet injections, as contrasted to injections with no momentum, 
can help to reduce the pipe flow distance required for the mixing of injected and 
 embiont fluids without requiring appurtenances or devices inside the pipe. For uniform, 
{ _ turbulent flow, experiments were conducted to determine the optimum injection 


: - conditions for a single jet at the pipe wall at 90° to 150° relative to the ambient flow, 


16539 EQUATION FORMULATION FOR GROUND- WATER FLOW 


and for two jets at 90°. Depending on the type of injection and the desired degree of 
mixing, the required flow distance can be reduced up to 70 percent compared to © 
injection with no momentum. For uniform flow, there are only small benefits from 
using more than two injection points. However, the two experiments which were 5 
conducted with a secondary current in the oe flow indicated greater benefits tl 


| REFERENCE: Fitzgerald, Steven D. (Sr. Grad. Engr., Turner Collie & Braden a 


: Houston, Tex.), and Holley, Edward R., “Jet Injections for Optimum Mixing in Pipe 
Flow,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY10, Proc. roc. Paper _ 


16535, October, 1981, pp.1179-11995 
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KEY WORDS: Boundary value problems; Boundary values; Computation; 
_ Computers; Ground water; Numerical analysis; Porous media; Porous media aietens | 


ABSTRACT: Advantages of | boundary tateaeit- equation methods (BIEM) for the 
analysis of two-dimensional flows through porous media are numerous. The data 
required and the number of unknowns are significantly fewer than for other general, 
numerical methods such as finite differences of finite element methods. A BIEM based | 
on the Cauchy integral equation is developed. This formulation has the advantage that 
= equation used for any unknown is a Fredholm equation of the second kind which 


— 


results in well-conditioned equations and allows nodes to be closely spaced in the 
region of singularities. The method is extended to the analysis of multizone, — 
anisotropic flows. Computational algorithms for solution of the equations by a simple 


i iterative method and by direct solution of the full set of simultaneous linear equations 


REFERENCE: : Hunt, Bruce (Reader i Civ. Engrg., Univ. of Canterbury, 
Christchurch, 1 New Zealand), and Isaacs, Lewis T., “Integral Equation Formulation — 
for Ground-Water Flow,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. 


HY 10, Proc. Paper 16539, October, 1981, pp. 1197- 
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ABSTRACT: An experimental study is described that concerned the surface erosion of _ 
an illitic silty clay (Grundite) at selected salinities and water contents; the study was a 
designed to test the applicability of rate process and double layer theories. Both — 
velocity increment and temperature increment tests were run in a refrigerated water 
tunnel. Experimental activation energies and flow volumes were computer from test 
results. Values for experimental activation energies range from 14 to 32 =s 
suggesting that interparticle contacts are essentially solid to solid as for soil 
deformation at higher stress levels. The mechanisms controlling resistance to surface - 
erosion are fundamentally similar to the mechanisms controlling soil — =e 
REFERENCE: Kelly, William E. (Assoc. Prof., Dept. of Civ. ioe Univ. of Rh Rhode 
Island, Kingston, R. I.), and Gularte,, “Erosion Resistance of Cohesive Soils,” Journal — 
_of the Hydraulics Division, ASCE, Vol. 107, No. HY10, Proc. Paper 16587, October, ‘ 


ie Equivalence; Flow charting; Flow rate; Free surfaces; Ice cover; Numerical “i oy 
ABSTRACT: The distributions of turbulent « eddy, and mass diffusivities were 
computed for three pairs of free-surface and ice-covered flows using the k-€ turbulence 1 
: £ model. The flow rate, channel slope, and bottom roughness values were kept the same 
| for each pair of flows and the flow depths were computed. The resulting flow depths oe 
: for ice-covered flows were 15 to 30 percent higher than the depths for free-surface _ 
flows. The computed velocity and eddy viscosity distributions do not follow the 
~ conventional logarithmic and parabolic distributions for the whole depth of flow. 
Concentration distributions resulting from the introduction of a & 
REFERENCE: Lau, Y. Lam (Head, Environmental Hydr. 
: ; Div., Nationai Water Research Inst., Canada Centre for Inland Waters, Burlington, 
; Ontario, Canada, L7R 4A6), and Krishnappan, Bommanna G., “Ice Cover Effects on 
| Stream Flows and Mixing,” Journal of the Hydraulics Division, a Vol. 107, No. 


| HY10, Proc. Paper 16602, October, 1981, pp. 1225-1242 i icap Aa 
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16603 COMPARISON OF FREE SURFACE VORTICES he 
KEY Field tests; Free surfaces; Inlets (waterways); Models; rth 
ABSTRACT: A Comparison between model and prototype observations of vortex 
| intensity is used to investigate potential scale effects of modeling free surface vortices. 
Due to the inherent difficulty with field observations, a reasonably large number of | 
such comparisons are needed prior to reaching conclusions. To this end, a survey of 
available model-prototype data is made and presented in summary form for a variety 
of intake types. Based on this comparison summary, and referring to previously — 
published research on this subject, , some scale effects may exist when modeling air core , 
vortices, and some increase in model flow may compensate for such effects. Little or 
no scale effect is evident for weak vortices or dimples. Recommendations are made for 
REFERENCE: Hecker, George E. (Dir., Alden Research Lab.; also, Assoc. Prof., 


Engrg. Dept., Worcester Polytechnic Inst., Holden, Mass. 01520), “Model- — 
Comparision of Free Surface Vortices,” * Journal of the Hydraulics Division, ASCE, 
107, No. HY10, Proc. Paper 16603, October, 1981, pp. 1243- 
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S. -Customary-SI Conversion Factors 


, _ SI (International System) units, the following list contains conversion factors t to enable readers ; 
nha to compute the SI unit values of measurements. A complete guide to the SI system and its y a : 
use has been published by the American Society for Testing and Materials. Copies of this a 
publication (ASTM E-380) can be purchased from ASCE at a price of $3.00 each; orders must A> 


a a” All authors of Journal papers are being asked to prepare their papers in this dual-unit om, 
>. provide preliminary assistance to authors, the following li list o of ene factors a and | guides * 


by the ASCE Committee on Metrication. 
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square inches (09 is in. 
square feet (sq ft) q 
square yards (sq yd) +; 
square miles (sq miles) i square kilometers (km* ) 

inches (cu in. cubic millimeters (mm 
cubic meters (m ) 
cubic meters (m’) 
kilograms (kg) 
kilograms (ke) 

kilogram force (kgf) newtons 

be per square foot pascals (Pa) 
pounds per square inch (psi) kilopascals 


acre-feet (acre- cubic meters (m 
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IN ‘CALM WATER 
‘By Jean L. Boillat® and \ and Walter H. Graf,’ M. 


During a study on the effect of turbulence on settling the 
writers (4), some interesting experimental observations have become evident 
on the. ‘settling velocity of spheres in a calm, thus nonturbulent, water. ‘The 


writers would like to communicate the these findings herewith herewith. 


liza. bas 4d) (5) 


The experimental installation was described in detail in Refs. 2 and 3; ‘it 


should suffice to give e but a brief review. . The spheres being , perfectly spherical 

_ and smooth, had diameters of 0.200 cm = D = 1.905 cm and densities -— 7 

.018 g/cm* = < 2.700 g/cm’. They were falling in a rectangular 
and transparent container having the dimensions of 130 cm x 130 cm x 80 > 
fd cm; the latter was filled with water of a temperature of 17.6° = T° = 24.1°, 
thus of quasiconstant density of p = 0.9983 g/cm’ and of a viscosity of 0.0106 
cm’/s < v = 0.0091 cm The experiments covered a Reynolds number 
range of 140 = R = 13,394, with R = v,,D/v. To ensure that the particle’s 
settling velocity, v,,, is not influenced by boundary effects such as the surrounding © 

7 & walls, the particle’s departure and arrival, a prestudy (see Ref. 2) was performed > 
which rendered that measurements should only be taken in a the zone 


Special attention was | taken ‘the ‘spheres w were re introduced, with 
_ departure velocity, without rotation and without air bubbles being attached. 
_ The spheres were released 5 cm or more underneath the ceiling, using a special 
Denier ity device to achieve this. In order to measure the time of settling 
oe. the following system was as developed: 2 cameras and a : video recorder 


Research Asst., Lab. d’hydraulique, Ecole Polytechnique Federale, (CH-1015 


? Prof. and Dir., Lab. @ hydraulique, Ecole Polytechnique Federale, CH-1015 Lausanne, 

_ Note.—Discussion open until March 1, 1982. To extend the closing date one month, 

a written request must be filed with the Manager of Technical and Professional Publications, 

_ ASCE. Manuscript was submitted for review for possible publication on December 10, 

4 1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of the 
American Society of Civil Engineers, @©ASCE, Vol. ‘107, be HY10, October, aww" 
ISSN 0044-796X /81 /0010-1123 / $01.00. 
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FIG. 1. —Schome of Container: “Showing DEC, Are Por 
in the measuring zone (DEC); and (2) the clock indicating entrance and a 
time respectively. The entire process was also restituated on a TV monitor. a 


"Subsequently the taped images were analysed by reading the time of passage 


of the sphere into and out of the measuring zone. In this w: way a precision 


the time measurement of 1/40 sec was achieved. 


the conventional drag 

‘coefficient, CT, versus Reynolds number, R, in Fig. 2(a). tess 
In Fig. 2, one immediately observes that the present data do- not fall on 
= “‘standard’’ CTS = f(R) curve as given, e.g., by Schlichting (Ref. 8, p. 
17). Since the experiments were performed under extreme experimental control 
the writers exclude, for the time being, that these deviations could be attributed 

_ However, it became evident that a | repeated experiment, with the same sphere 

_ and under identical experimental conditions, showed a certain scattering. It 
is therefore that the ‘‘median’’ of each set is more representative; see Fig. 


2(b). In Table 1 these ‘ ‘median’ data are with all the experimental 


A polynomial regression fit to the present data toa new ‘*reference”’ 
CTR = JS(R) curve, which is given in Fig. 3. The dotted lines, at the CTR i 
_= f(R) curve on the two extreme ends of the experimental range, indicate 

_ that due to lack of sufficient data, a certain caution is in place. Nevertheless — 

_ within the zone of 2 x 10’ < R < 6 x 10’, the ratio of CTR/CTS is always 


higher | than 1.10 and may get as high as 1. 3. Furthermore, it will not escape _ 
to the observer that in the CTR plot at R < 2 x 10° there exists 


—_ 


_ Before the writers’ results are considered, it seems in place to review whether 


i _ similar observations on deviation from the standard curve have been made t before. I 
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free settling ina quiescent fluid, 
ly The research available before 1930 is well-reviewed by Schiller (7) and is 
= _ Teproduced in Fig. 4. There are the e freely-falling spheres in laboratory experiments 7 
* by Allen (1900), Schmidt (1920), . Liebster (1927) and Lunnon (1928); there are 
is 4 the freely-falling spheres in outdoor conditions by Shakespeare (1914), Lunnon 3 
= (1924) and Bacon and Reid (1924); there is the freely-rising balloon experiment © 
_ by Hesselberg and Birkeland (1912) and there is the sphere’s drag experiment . 
in a windtunnel by Wieselsberger (1914). In addition it should be said om 
Schiller (7) was careful in choosing only data of high experimental quality. 
shows that above R deviations though not always 


literature, there are only few investigations available which in fact study | 


RE 


_ systematic from the standard CTS = = f(R) do occur; even a certain hump is_ 
“ In a study reported by Stringham, et al. (9) it is also evident that the “fixed — 
sphere”’ data and the “‘falling sphere’’ - data give different results. This may 
be seen in Fig. 5, where the present Laboratoire d’hydraulique, Ecole ~~ 
nique Fédérale (LHYDREP) data are plotted as well. Interestingly enough the el 
hump appears to coincide with our experiments. _ Stringham’ s et al. (9) data 
“were obtained with spheres of steel and plastic, which fell freely in different 

_ liquids such as water, water-glycerine mixtures and pure glycerine through 4 
cylindrical tube, 40 cm in diam and 300 cm in height. These authors observed 
_ that notably the plastic spheres “‘exhibited small erratic horizontal movement.’ 
the authors assign the Govtetions of their dota ‘from 


a 
FIG. 4.—Data from Schiller (7) Compared with Present LHYDREP Data 
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_ The present writers submit to have no reason that their experimental setup — 


~ and the data- -taking procedure can be made responsible for such large deviations — 


i "deviations from the standard CTS = = f(R) relationship are after all also ‘reported me 
Schiller (7) and more recently by Stringham, et al. (9), 
_ In the following the writers shall put forward some speculative ideas as to 
_ : why the present data’s CTR = f(R) show differencies if compared with “‘standard”” 
- Firstly it appears that the standard CTS = ff (R) curve was established mainly—if_ 7 
‘not exclusively—with data obtained on fixed spheres i in wind tunne els; and this 
: _ particularly i in the region of R = 800 [see Murray (6) p. 304]. Already Schiller 
ma) pointed to the fact that a windtunnel has always a certain degree of turbulence 
a imposed upon the free-stream velocity, which could modify the drag 


coefficient. Secondly, it is important to note that in windtunnel experiments 


FIG. 6.—Data from Second (5) for Spheres and Disks Compared with Present 
‘the spheres are fixed in baie However, a freely-falling sphere will adjus 
_ its path according to R dependent vortex shedding at its rear, which in tun 
will lead to an oscillatory descent of the sphere. Theory of similarity a tn 
suggest that an adimensional frequency of oscillation—such as the Strouhal 
number S = nD/v, ,—Should become an additional parameter. Unfortunately 
the measuring technique was not designed for observation of the | entire ‘path 
of the falling sphere and thus for evaluation of these frequencies. On the other 
+ a good S = f(R) relationship for spheres is at the present not available _ 


[see Achenbach (1) and Taneda(10)}). 


From the foregoing it is reasonable to expect that spheres falling freely through 
calm tender different results. Interestingly enough and thus 


4 


was - obtained for freely-falling disks and ‘does show : a 4 pronounced hump. A 

falling disk wabbles, resulting i in a pronounced Strouhal number effect. While 

there exist to the writers’ knowledge no fixed- disk data in this hump region oa 


_ 
— 
4 &§ 


1 

(R = 300), Schiller (7) an of a a would 


The settling velocities « of round and smooth spherical particles in calm water i= 


10° . < R< 6 x 10° were analyzed utilizing modern measuring techniques. 
If the present results are compared with the “standard” CTS = f(R) curve, 
_ they fall on the average 20% higher as is to be seen in Fig. 3. However, the — 
_“‘standard”’ curve was established with the majority of data derived from fixed = ; 
- spheres tested in windtunnels. The experiments of Schiller (7) and of Stringham, 
et al. (9) show the same tendency as ; observed in the present data. To make » 3 _ 
2 ‘a recommendation for a new drag- -coefficient curve for freely falling spheres . 
is premature in the light of the writers’ limited experiments. To settle the existing — - 
discrepancy between “*fixed”’ and “‘falling’’ spheres, it will be necessary that 
further research is dene, possibly focusing on the effect of the ‘Strouhal number 
8 on the drag of of freely falling spheres. can 7 4 
my 
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; The following symbols are used in this paper: 


CT = "general drag coefficient; 
CTR = reference drag coefficient; Sait 
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SPHERICAL PARTICLES 


= diameter of spherical pesticle; 
nD/v,, = Strouhal number; 
v,, = spherical particle’s settling velocity; 
of water; and 
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LENGTH OF FLow SEPARATION OVER | Dunes 


of the interrelationship between the fluid and the bed. When the flow is tranquil 

; a (Froude number <1), the water may deform the bed into ripples or dunes depending ¢ 

1 on the flow conditions and the grain 1 size of the bed material. As the = 

_ passes over each bed form, separation occurs at the crest creating a wake mn 
or eddy which extends some distance downstream to the point of reattachment _ 
on the back of the next bed form. Within these separation zones, reverse — 
circulations exist. Knowledge of the length of the flow separation is useful — 
in solving a variety of f hydraulic problems. Crickmore (2) and the writer pe 
Lau (3) developed equations to compute bed- -load “discharge from bed -profile 
records of migrating dunes, in which the position of the point of reattachment — 
was an important consideration. Vanoni and Hwang (16) developed a relationship b 
_ for the friction factor in which the projection of the upstream face of the 
bed forms exposed to the flow must be known. Chang (1) also « derived a friction 7 
_ factor equation, valid for r ripples, which required the length of the eddy | ena ‘ 
: behind the bed-form crests as one of the variables. Pillai (1 1) presented a correction © 
” factor for the average flow depth over a bed composed of sand waves, the 
se of which varied with the steepness of the bed forms. It is most 


§ probable that this depth adjustment is related to the length of the flow separation, © 


thereby defining the average depth of the water mass which is moving downstream. 
~The e writer and Lau (4) and © Vittal et al. (17) found that it is. not correct to S 
= that the friction factor due to the sand-grain roughness on the back 
of a dune is the same as that due to a plane bed of similar grain size. This 
| probably also due to the fact that the length of the bed form exposed to 
the flow varies with the dune and ‘must therefore also be ‘Telated 
_ to the length of the flow separation. 


‘In this paper, an investigation of the _e of the “flow, panels behind 
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dunes is made. Consideration is given to 0 the effects of the dune geometry, — 
-sand- ~grain and defining the flow conditions. 


oar review of the literature showed that attempts have been made to obtain 
& length of the flow separation for different types of bed topographies. Raudkivi — 
& 14), in studying flow over ripples, found that in the lee of a ripple there 
_ is a much stronger agitation of the sand surface at five ripple heights-eight - 
( heights downstream from the crest. Walker (18) measured the separation — 
length behind a rearward facing rectangular step and found it to be slightly — 
i scccups) than about five times the step height. Tani (15), using a similar step, 
- found the eaperatien length to be about six times the step height. Measurements 
_ by Etheridge and | Kemp (6), also over a rearward facing step, showed that 
“the ratio of the separation length to ‘step height varied between 5 and % “Chang. 
(d+) measured the length of the separation zone behind a single cast iron angle z 
placed on a plane bed in a rectangular flume, and repeated his tests for three 
different sizes of angles and three different flow conditions. He found that 
a the ratio of the separation length to the height of the angle varied with the 
- flow velocity from 5- 8. ‘However, it is not likely that flows « over a a single Lad 


Knoroz (10), i in developing an n equation for the friction factor for dune beds, 
4 determined that the length of the separation zone is ten times the dune height. 
Examination of data by Jonys (7,8) indicated that for his dunes, the length 
the separation zone on the average about 4. 2 times the dune 


the ratio of separation length to dune height w was constant and ‘equal to » 3 3, 
They used a new technique which utilizes the streaming birefringence of dilute _ 
colloidal suspensions of birefringent materials, 

_ The results from the literature indicate that the length of the separation zone 
be behind dunes may vary from five times-ten times the dune height. Unfortunately, 
the information is too limited to determine a general relationship between the 
; separation length, bed- form shape, sand- -grain roughness, and flow conditions. - 
_ The results from measurements over dunes and ripples and rearward facing 
‘ ‘steps appear to be similar. However, it is doubtful that this similarity can be 
taken to be general, especially when the sand waves have steep upstream faces. 
In view of the limited information on the flow separation length, a set of 
experiments was conducted in an attempt to obtain a comprehensive relationship 
_ between the length of flow separation over d dunes and the pertinent independent 

. In order to have good control of the variables defining the dune height and 

~ length, it is convenient to use a rigid bed composed of artifical dunes. Under 
these conditions, considering two-dimensional flow in a straight uniform channel, . 
the flow separation length can be to Fig. l(a) by 


a set set of independent veer 


| | 


height © of the bed forms; @ = angle of inclination of downstream face of the 
< dunes; D,, = median grain-size diameter; h = average depth of flow; coll 
\ Sem flow velocity; p = density of the fluid; » = viscosity of the fluid; 
S and g = acceleration due to gravity. Eq. | om be written in dimensionless 


, for >0.. mip 


Measurements by Etheridge and Kemp (6) have: shown that flow in the 
‘separation zone is quite insensitive to the flow in the corner region at the 
foot of a reverse rectangular step. This suggests that quite crude — 
of the dunes in this region between the | crest and downstream toe is justifiable Bh 
Therefore, since ‘the angle does not normally "vary a great deal, it is not 


a significant pear in Eq. 28 = can be omitted from further consideration. 


J FIG. 1.—(a) Definition Sketch for Dunes; | (b) Definition Sketch for Reverse Steps -_ ; 


n the dune regime, of the Reynolds. number are so that 
: o effect of viscosity is also not significant and the Reynolds number given 


Uhp/p can be deleted. The relationship being ‘Sought i is. then wy 
= —,- j for A 
of > 


Experiments: were to. investigate the variations of L/A 
independent variables in Eq. 3. vornd mow diced al owl 


: Experiments v were conducted in a tilting flume rectangular in cross section, — 

20 m long and | m wide. Discharges were measured at the downstream end * 


using ® weir box. The flow leaving the flume was split into three parts = 


| 
ont 


— 
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“portion of ‘the | flume, the effect of the ‘side walle ‘was practically eliminated. 
The flows were always uniform with Froude numbers not exceeding weeraerers 


_ The dune roughness was obtained by using artificial bed forms fabricated — 
from plywood and covered with sand. Four different bed-form shapes, — 


sand was fixed to the bed forms by first covering these with “Mactac,”” 
self-adhesive vinyl covering. Then a thin layer of varnish was ‘brushed on, 
~~ and the sand was sprinkled on | Lop so that the ‘ “Mactac” was completely —— 
(11). After: drying overnight, the ¢ excess sand was washed off and a 


= 108 pos of & = = 30°, and were e placed aby the full length of the flume. ‘iT he 


the “Mactac””_ was s stripped off and the process ‘repeated. Each bed form v was cy 
tested with one or more sand sizes before a new one was placed in the flume. 
In all, three different sand sizes were used having median diameters of 0.62 — 
7 _ ‘mm, 1.20 mm, and 2.60 mm. Therefore, each sand size was approximately 
twice as large a as the. one preceding it. Tests on the smooth bed forms were — 
eonducted using the smooth surface of the ‘‘Mactac’’ ‘only. A special set of 
tests was conducted in which a rearward facing step with horizontal upstream 
and downstream boundaries was used, as shown in Fig. 1(b). The height of | 
the step was the same as that of the dunes. This configuration was tested a 
with a smooth surface (i.e. a o ™ 0), when covered with a sand of 1.2-mm _ 
P median diam. - This set of tests was done to provide | data for comparing the 
_ length of separation zones behind the crest of dunes and such a reversed step. 
m After examining the methods of measuring the separation length used by 
Chang (1), Etheridge and Kemp (6) and Karahan and Pederson (9), it was ~a 
decided that a simple flow visualization method would be adequate for the 
} present study. The flow separations in all cases were measured using a Potassium — 


Permanganate solution as a visual aid. A 100-ml hypodermic syringe with a_ 
_ needle 45 cm long and |-mm inside diam was attached to a traversing mechanism _ 
- on a movable carriage which traveled on the rails on top of the flume “ie. 
_ This permitted movement of the syringe in both the vertical and horizontal — “q 
_ direction. The tip of the hypodermic needle was ground off at an angle of 
—-30° with the » vertical. This made it possible to place the needle so | that ‘the 
ee was always injected horizontally at right angles to the main flow, and 
ensured that no momentum in the streamwise direction was imparted to the 
dye. Once a given flow was set up, the syringe was filled with the dye and 
positioned on the center line of the flow, and then lowered so that the needle 
was about | mm or less above the dune surface. 
& ‘The location of the point of flow reattachment was obtained by the ne following a 
~ two steps. In both steps, dyes were injected into the flow in small short pulses, = 
and were observed to note their direction of movement. In the first step, the 
Syringe was always positioned well downstream of the anticipated —— 
zone, so that there was a clear and distinct movement of dye pulses downstream 
_ [Fig. a)]. The syringe was then moved in small increments to new positions 
upstream, and in each case, the behavior of the dye pulses was carefully noted. 
‘The progressive upstream positioning continued until it was observed that some 


| 


of the dye tended to trav 
edge of the region in which the | point of flow reattachment occurs. The concept = 
of a region was adopted since, due to the turbulent fluctuations (6), it would | 


not be possible to define a precise point, but rather an ae In the second - 


the separation zone being considered. At this point, the dye Pao te clearly moved — 
ap upstream and up the s steep face of the dune. The syringe was then placed — 
at different positions, each being about 3 mm-4 mm further downstream than 
the previous one, and the movement of the dye observed [see Fig. 2(b)]. This 
was continued until once again a portion of the dye pulse was entrained in 
the downstream flow. This point was taken to be the upstream edge of = 
region for the point of flow reattachment. All positional measurements were 
pr to the crest of the dune upstream of the separation zone and thus, 
the probable location of the of flow was 


énebde 


Point 
_ the average of the measured pairs. Several pairs of aiads measurements were 
made for a given flow condition to ensure a reasonably reliable position for 
the point of reattachment. The variation between values of the length of flow 
; z separation» for the different pairs of measurement was less than 15%. _ These 
“results are are encouraging when compared with the results of Etheridge « and Kemp 
- (6). In all, data for 291 runs were obtained (4), in which dune length, flow 
a depth, sand size, and the discharge were the basic variables that were varied. 


ad Separation Length behind Dunes.—Values of L/A were plotted versus the 
Froude number, F = U/V gh, for each dune shape, using h/A as 
_ parameter for each given value of D,,/A. Due to the large amount of date, 
_ only the plots for A/A = 0.02 and 0.07 are shown as Figs. 3 and 4. The 


computed values of F, are om in Ref. 5. In all cases, 
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“of the Froude number, and is, at best, only weakly dependent on h/A. "Therefore, 
since F and h/A are shown to be insignificant, they can now be omitted from 
further consideration. The data analysis is then reduced to searching for the a 
A toe cate for the dunes, aod thal of 


Average values of L/A were obtained by fitting smooth straight lines through 
bes: on the plots of L/A versus F= = U/V gh f for all values of A/A and Dz,/ : 

= ested. Typical examples of this ; are also given in Figs. 3 and 4. The e values 
of L/A thus obtained, together with the corresponding v values of A4/A and © 


— Dso/ /4 A, are given in Table 1. The values of L/A were then plotted as a function — 


FIG. 4.—Variation of D Dimensionless Length of Flow Froude 


When A/A=0.07 
ia 


of A/A using D,,/A as a "parameter in Fig. 5. The relationship in Eq. 4 is 
_ Shown as a family of diverging nonlinear curves with one curve for each value | 
D/A. For each curve, values of L/& decrease as A/A increases from 
_A/A = 0.02-0.07. The rate of change of L/A with A/A appears to be greatest 7 
when A/A = 0.02, decreases progressively as A/A increases, and becomes _ 
very small when A/A = 0.06. The curves also show that for small ~ ll 
of A/A, the effect of change in grain size is more significant than for larger — 
of A/A, and when A/A = 0.06 the effect of D,, becomes virtually 
should be noted at this” point that the test idealized in 
_ Many respects in order to provide a wide enough range of data to assess Ne 4 
; _ effect of the important independent variables. In reality, values of A/A in the 
will probably not t exceed 0. 06 (19), and values of will 
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ratio L/A was again independent of U/V gh and h/A. When pion the 
Brg step, it is useful to regard it in terms of the dune shape parameter — 


TABLE 1.—Average Values of te 
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@ Raudkivi (13, 14) 
Jonys (6, 
@,4,8 Present Data 


A - > and for this case RMA =0. 


Doo 
Eq. 5 shows that for the reverse step, L/A is a function only of the surface 
¥ roughness of the step, a result which has not been reported in the literature. y 
_ This may partly explain the difference in reported values of L/A for such 


= 
a Separation Length behind Reverse Step.—The data for the reverse step were 
reated in ne ame ne Cine state ne ate On ath ne 
| 
“a 


Values of L/A we were plotted for A/A = 0, and the two cases of D,./& = ~& 

0 and D,,/A = 0.040 in Fig. 5. The values of L/A for each case are higher 

than the values for dunes for the same values of D.,/ A, but the trend of decreasing 

L/4A with increasing D,,/A is still maintained. The curves for values of D,,/A — 
= 0 and D.,/A = 0.040 were smoothly connected to the points of L/4 for 

the corresponding. values Of Dso/ A at A/A (see dashed segments in Fig. 
5). The curves from A/A = 0-0.07 are interpreted to mean that the vertical — 
step is nothing but the limiting case for the dunes, and that values of L/a 
7 for dunes (A/A > 0) can never be greater or equal to values of L/A obtained 
fo a reverse step for a given value of D,,/A. Therefore, any comparison of — 

L/A between du dunes and reversing steps is incorrect, and errors incurred incremse : 


Comparison of Results with Data from Other Sources. —Karahan and Peterson 
(9) obtained a value of L/A = 5.3 for their dune shape of A/A = 0.06. This 
_ value was plotted on Fig. 5 which shows that their value of L/A is considerably _ 
higher than the corresponding value obtained from the present tests. Karahan — 
and Peterson (9) stated that their value of L/A agreed well with the a 
obtained by Tani (15) and Walker (18) using rectangular reversing , steps, as 
well as those results obtained by Raudkivi (13) using a fixed ripple bed form. 
It was shown that comparisons between such reversing steps and dunes are 
invalid and in the case when A/A = 0.06, values of L/A for reversing steps — 
_ are much higher. Therefore, the favorable comparison with the rectangular step - 


On the other hand . the ‘ripple form used by Raudkivi (13) 5 was 5 steeper ‘than 
A/A= 0.06, having a value of A/ A = 0.075. Karahan and Peterson (9) computed 
a value of L/A = 0.32 from Raudkivi’s data which would —_— in a value 
of L/A = 4.3 for A/A = 0.075 (i.e., L/A = L/A-A/A). This value of 
L/h was" also plotted on | Fig. 5, showing good agreement with the he curves pm 

_ Jonys (7) conducted a series of tests over mobile beds, for which the average 
valees of A/A and D,,/A were about 0.06 and 0.06, respectively. The position 
of flow reattachment could be obtained from -Jonys (8) and the value of L/A — 
‘was computed to be » 4.2. This value of L/A was plotted on Fig. 5 at  A/A 
= 0.06. Considering that D,,/A ~ 0.06, the point should plot between the curves 4 
for D,,/A = 0.040 and 0.087. However, considering the variance in mobile | 


- boundary measurements, the agreement of L/A from Jonys (7,8) with that from 


& the > present tests is quite good. 


‘The ony of the flow separation on the lee side of dunes of dif ferent echaliliiien’ 
_ dm roughness, as well as for a reverse rectangular step was measured and» 
the following conclusions were reached. 
-_ The length of the flow separation behind dunes is independent of the Froude 
‘aember and virtually independent of the flow depth for Froude numbers a 
than 0.5 and values = 
2. The length of the flow separation behind dunes is dependent on the sand- grain j 
"roughness on the dunes when the dunes are flat. ‘When dunes are (i. €., 


| 


MMe wor 


A/A > 0. 05), the sand-grain ‘roughness effect is negligible. In fact, 


the separation length for steep dunes is approximately equal to four times the > 
3. The reverse rectangular step may be considered to be a special case of 
dunes with infinite wave length (i.e., A/A = 0). The largest attainable flow 
a "separation occurs for this case. It is therefore not correct to compare flow 
separation lengths behind dunes (ie Oui > 0) directly with those behind 


elec 


_ 4. The length of flow separation behind dunes for a given surface roughness | 
increases significantly for values of A/ A < 0.05 as A/A decreases. This variation -_ 

most significant for values of D,,/A ~ 0. As D,,/A increases, Tate of 


change i in length of flow separation becomes less as A/A decreases. _ feast DVB = 
5. The length of flow separation behind a rectangular reverse step is a function a 
surface roughness only for the range of flows tested. 


6. Fig. 5 can be used to obtain estimates of the length of flow ‘separation: 
for most situations with dunes. bes nyt ry 
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The following are used in this paper: 


Doo = median grain- size diameter; 
4 g acceleration due to 
= average depth of flow; 
length of flow separation; 
average velocity of flow; 
length of dunes; 
p= viscosity of fluid; 


angle of of downstream face of dunes. 
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Don J. Wood,' M. ASCE and A. G. 


Steady- state of pressure and flow ia systems is problem 


a great importance in engineering. The basic hydraulics equations describing — 

: the phenomena are nonlinear algebraic equations which cannot be solved directly. 

_ These equations have been expressed in two principal fashions. They have 
been written in terms of the unknown flowrates in the pipes herein referred — 
to as loop equations. . Alternately, they have been expressed i in terms. of unknown 
“heads at junctions throughout the pipe system (n (node » equations). Several algorithms * 
have been proposed for solving these equations and these techniques are _ 


_ In this paper the principal algorithms are presented for applications to pod 

systems of general configuration and systems which contains pumps and other 

_ common | hydraulic components. The reliability of these algorithms a are : examined — 


; by solving a large number of actual pipe network pr problems \ with each ) algorithm 
_ There is a considerable amount of published material « dealing v with pipe network 
= all of which cannot be reviewed here. However, some ol 
contributions of historical interest will be cited. Hardy Cross authored the original — 

and classic paper (4). In that article, which considered only closed loop networks — 


at no pumps, a method for solving the at equations based on adjusting : 


_ Cross method. Although it is not as widely used, Hardy Cross described a 
second method for solving the node equations by adjusting: the head at each 
node so that the « continuity equation is balanced. A number of subsequent p papers 
have appeared which further describe these methods or computer programs ~ 
utilizing these methods (2,5,8,11,15). 
_ Because the adjustments are computed independent from each other, conver- 
gence problems using the methods described by Hardy Cross were frequently 
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i _ This method is very widely used today and is often referred to as the Hardy “a 
| 
| 
‘4 


noted and procedures developed to improve convergence. Martin and Peters vebl ' 

_ (14) and Epp and Fowler (6) described a procedure to simultaneously compute aes / 

_ the flow adjustments for closed loop systems. This procedure had much improved — 

_ convergence ¢ characteristics and forms the basis for more general applications — 
‘Both the methods just described for solving the loop equations require an : 

initial balanced set of flowrates and the convergence depends to a degree on 

how close this initial set of flowrates is to the correct solution. A third method, = 
developed by this writer, solves the entire set of hydraulics equations simulta- | 


neously after linearizing the nonlinear terms in the energy equation. This procedure 


was first described for closed loop systems (18) and has subsequently been 7 
modified for more general applications (7,17). A similar approach has been 
developed for the node equations where all the node equations are linearized 
= solved simultaneously, and this method is described for closed loop systems — 
by Shamir and Howard (16). ‘Abtiiaeal references to this method have been 
"The five methods just pang described in the next section for 
applications. No restrictions on network geometry are required and pumps may > 
be included anywhere in the network. Certain special components such as check © 
valves and pressure regulating valves are not. considered since these components 
require special treatment which depend or on algorithm employed. Each of 
these methods require iterative computations where the solution is improved 
until a specified convergence criterion is met. If a sufficiently stringent conver- 
gence criterion is met, the solutions normally will be essentially identical for 
all five methods. In some cases, however, it is not possible with certain of 
= “the algorithms to meet a specified convergence criterion ‘regardless of thet number : 
_ of trials completed. In other cases a seemingly stringent convergence ce criterion 
may be met but the solution is still in considerable error. Convergence difficulties 


as these have been previously noted and reported. 


In the original paper by Hardy Cross (4) he noted that “‘convergence was 
3 slow and not very satisfactory” when employing the head adjustment method 
he. developed. This was attributed to ) using initial head estimates which were 
“not very good. Of the two methods “described by Hardy Cross, the method 
of adjusting flows became the most widely taught and used method. Convergence — 
een using this method were also recognized, however, and several sugges- 
tions were made for i improving convergence. Investigators have advocated the 


and “suggested using | a selective procedure for choosing loops. 
= as a means of accelerating convergence (10). It appears that these and other 
procedures suggested for improving convergence of this method will improve 
only in certain situations and will not assure convergence. 
Convergence problems are largely unreported for the improved methods 
developed for solving the loop equations. However, additional convergence 
problems have been reported for methods based on the node equations since | = 
Hardy Cross first alluded to such problems in his original paper. Dillingham -_ 
_ (5) stated that when the method of adjusting heads at individual nodes is applied — 
to a large network it may not converge or may converge very slowly. He described | 
- some procedures for improving convergence. Robinson and Rossum (15) who | 
developed a computer program based on this method state that, “‘convergence ~ 
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They further state 
that ‘‘convergence may not occur if check valves are present.” method 
of simultaneously adjusting heads normally converges much more aagidiy which > 
lead Lemieux to state that convergence was assured (12). However, it appears 
that this assessment is optimistic and problems have been noted with this me method. 

= have been noted by Shamir and Howard (16) who also report that 
there is a possibility that a solution can not be obtained. Liu also stated “for 


poor initial input the method (simultaneous node adjustment) may diverge from 


the true solution or converge slowly’’ (12). Collins and Kennington presented 
some data which documented _ convergence problems f for a large network using 
this method (3). ‘the unkdown he ade ai junetion a 
The reliability of the algorithm employed for pipe “network | analysis is hie i 
great importance. Failure to obtain a solution is an inconvenience and the feiiess 
to recognize a poor solution may be even a greater problem because this may 
lead to poor design or management of water distribution systems. The i 


Atcorrrums For Pipe Network ANALYSIS tS 
The commonly used algorithms for pipe network analysis are described ine 
this section. Detailed presentations can be founded in the references cited for 

4 each method. A brief description of pipe network geometry and the basic hydraulic 
4 


ae are included so the algorithms os be presented using a common 


‘tw Pipe Network Geometry.—Basic geometry considerations for a pipe network 

are summarized as s follows. A pipe network is s comprised of a number of pipe 

sections which are constant diameter ‘sections which can contain p pumps and = 
fittings such as bends and valves. The end points of the pipe sections are — 
y nodes which are identified as either junction nodes or fixed grade nodes. A © 

yj node is a point where two or more pipe sections join and is also 

a point where flow can enter or leave the system. A fixed grade node is a 
point where a constant | head is maintained such as a connection: to a 
tank o or reservoir or to a constant pressure region. In addition primary —— 

can be identified in a pipe network and these include all closed pipe circuits _ 
fd within the network which have no additional closed pipe circuits within them. a 


in which p = number of pipes; j = number of junction nodes; l number 


primary loops; /= = number of fixed grade nodes. 
a Basic Equations. —Pipe network equations for steady- -state analysis have been 


commonly expressed | in two ways. Equations which ‘express mass continuity ; 


and energy conservation in terms of the discharge in each pipe section have 
- been referred to as loop equations and this terminology will be followed here. 
A second formulation which expresses mass continuity in terms of heads at 
junction nodes s produces a a set of equations referred to as node equations denn 


_ is slow when a network cont 
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offers a basis for a set of 


equations to describe a pipe system. In terms of the: unknown discharge in 
each pipe, number of mass continuity and energy conservation equations 
a can be written equaling the number of pipes in the system. For each junction _ 
> node a continuity relationship equating the flow into the junction (Q,,) to the 


og Here Q. represents the ¢: external inflow or demand at the junction ‘node. For 
each primary loop the energy conservation equation can be written for pipe — 
sections in the loop as follows: shoo 


zh, = 2E, (/ Jequations) . 

in wee: h, = energy loss in each pipe section including minor losses; es 

E, = energy put into the liquid by pumps. ps. If there are f fixed grade nodes, 

Fa — | independent energy conservation equations can be written for paths of 

pipe sections between any two fixed grade nodes as follows: 


AE = 2h rh, - (f-1 equations) (4) 


in which AE = the head difference between the fixed grade nodes. The terms dd 
in the energy equations can be expressed as function of flowrate. Thus, the r 

_ continuity equations (Eq. 2) and the energy equations (Eqs. 3 and 4) form 

; ‘the set of p simultaneous equations in terms of unknown flowrates which are 

_ termed the loop equations. Since these are nonlinear algebraic relationships, _ 
fo direct solution is possible. Three algorithms for solving the loop venoee atl 
“are presented in this paper. Ke af the 
_ Node Equations. —The analysis is carried out in terms of an unknown total 
head, H, at each junction node in the piping system. The basic relationship _ 
used is the continuity relationship (Eq. 2). The flowrate in a pipe section connecting d 
nodes labeled a and 5 is expressed in terms of the head at junction node a, 

ay the: head at the of the pipe section, ane loss 


expression assumes that the pipe section contains no pumps and a head 


lationship is used of the form 


in which the term K = the loss coefficient for the pipe section and is a function 
of pipe parameters and flow conditions and depends o: on the head loss expression 


expresses Continuity at junction node 4, in which pipes connect, 


head loss expression 1 used. Combining Eq: Eqs. 2 and 56 gives: 
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of the adits in the summation depends on whether the flow is into or or out of 
> junction node a. A total of j equations are written in this manner. Gitiw Wrdreet 
_ The basic set of equations can be expanded to incorporate pumps. For each 
_ pump, junction nodes are identified at the suction and discharge sides of the 
-_ _ pump. Two additional equations can be written as terms of the two additional — 
~ unknown heads at the suction and discharge sides of the pump and the adjacent 
heads. One equation « expresses flow continuity in the suction and discharge 
lines (using Eq. 5) and a second equation relates the head change across the 
o pump to the flow in either the discharge or suction line. 2s” 
For a pipe network eds junction nodes and nm, pumps a set of j + 2n, 


which junction nodes and 
pump suction and discharge | heads @ at t all ‘pumps in the f pipe system. Like the 
_ loop equations, these are nonlinear algebraic equations and no direct solution — 
Algorithms for Solution of Loop Equations.—Three methods for the solution ‘: 
of the loop equations have been developed and are in significant use today. 7 
‘ Each use gradient methods to climinate nonlinear 


Adjustment (PA TH) Method. —This method of station was 
described by Hardy Cross (4) and is the oldest and most widely used — 
“The original 1 method ¥ was, }» however, limited to closed loop systems and included — 2 
a solution is summarized as follows. 


a Determine an initial set of flowrates which satisfy continuity at at each junction | 4 


= a “Compute a flow adjustment factor for the path of pipes ; for each energy yO 
equation which tends to satisfy the energy equation written for that path. The 
i; application of this correction factor will not disturb ¢ the continuity balance. — 
—- Using improved solutions for each trial repeat step 2 until the average 
correction factor is within a specified limit. T 
| The flow adjustment factor for a path i is computed from t the er energy iin 5 
_ for that path and is intended to correct the flowrates so that the energy ogres 
is satisfied. However, a correction for a particular path will disturb @ ‘enerey 


this method requires a flow adjustment to all paths pipe network (lo 


1 paths between fixed grade nodes). a pin 

Simultaneous Path Adjustment (S-PATH) Method. tn “order to improve 

; a convergence characteristics a method of solution was devised which — 
4 ~ neously adjusts the flowrate in each path of pipes representing an energy rey equation 
(6). This n method can | be be summarized as 
bagi Determine ar an initial set of flowrates which satisfy y continuity at each junction 

‘Gant compute a flow adjustment factor for each path which | 


- tends to satisfy the energy equations without disturbing the continuity balance. ie 


| 


Using the improved until the average flow 
The simultaneous determination of the path flow adjustment factors requires 
the simultaneous solution of / + f — 1 equations. Each equation accounts for 
the unbalance in the energy equation due to incorrect values of flowrate and 
includes the contribution for a particular path as well as contributions from 
all other paths which have pipes common to both paths. 
_ A set of | + f — 1 simultaneous linear equations are formed in terms of 
“flow adjustment factors. These linear equations can be solved using standard 
qi ‘procedures and the solution provides an improved set of balanced flowrates 
. | which can be used for another trial. Trials are repeated until a specified ac accuracy 
Linear (LINEAR) Method.—This method is based on a a simultaneous solution : 
of the basic hydraulics equations for the pipe system and has been reported > 
for closed loop systems (18) and general systems (17). Since the energy equations 
_ for the paths are nonlinear, these equations are first linearized in terms of 


e: equations (Eqs. 2) to form a set of p simultaneous linear equations in terms | 
_ The technique used to solve equations follows. Based on an arbitrary initial 

; value for the flow in each line the linearized equations i are solved using routine — 


to linearize the energy equations and a a second solution i is obtained. The procedure 
is repeated until the change in flowrates obtained in successive trials is insignifi- 
_ Algorithms for Solution of Node Equations.—Two methods for solving the © a 
node equations are widely used and these are described here. in 
Single Node Adjustment (NODE) Method. —This method was also described 
in the paper written by Hardy Cross (4). However, the method has never been 
as widely used as the single path adjustment method. It is, nevertheless, in 


l. prares a sian head | for each junction node in the system. . This g 
assumed head does not have to ) satisfy any conditions. However, the call 
© the initial assumption the fewer the requiredtrials, 
_ 2. Compute a head adjustment factor for each junction node which tends 
3. Repeat step 2 until the average correction factor for heads is within a > 


specified accuracy or some other specified criterion is satisfied. 


pe 


. = The head adjustment factor i is is the change i in head at a particular node which ad 
will result in satisfying continuity (Eq. 2) considering the heads at the adjacent 
nodes as fixed. Again a gradient approximation is used to compute the required salad 
; head change. A single trial for this method requires the adjustment of the head - 
for each junction node within the pipe system. The trials continue until the 
specified convergence criterion is met. 


Simultaneous Node (@NODE) Method.—This method is based 


4 
} 
J 
Signilicant use today. ihe method is summarized as = 
1 
| 
= 


ona a of the basic p pipe network and requires 

a linearization of these equations in terms of approximate values of the head 

: 4 oer This produces a set of j + 2m, simultaneous linear equations (in which 
, is the number of pumps). These “equations are solved as follows: sire 


set of | junction n node heads. ‘These 
- heads ar are » then used to linearize the set of node equations and the procedure 
repeated until Subsequent calculations satisfy a stated convergence or ecouracy 


A comprehensive comparison of the algorithms was an 
extensive data base was available. These data were provided mostly by consulting 7 
and water distribution engineers and represents actual or proposed 
distribution systems from all over the United States. In some cases these data 


sent to the writer because analysis difficulties were e encountered. These 
_ systems include very high or low resistance lines and pumps operating on steep 
curves. In general, however, the data is representative of the type of systems 
which are routinely analyzed. The data utilized includes 30 systems with under 
100 pipes and 21 systems with over 100 pipes (up to 509 pipes). The data 


for a _ System | can include changes which ‘Fesult i in the analysis: of additional 


into . the data. A total of 60 situations were analyzed for systems of under — 
100 pipes and 31 situations for systems of over 100 pipes. _ aismixorgge bis x 
3 _ The data includes pumps described in two ways. Pumps operating at constant — 
power and pumps operating on a curve passed directly through three points _ 
g operating data (head © vs discharge) were considered. Components such as 
of op valves and pressure regulating valves were yere not included since special — 
q procedures are required to handle these components. The Hazen Williams equation ; 
, for head loss was employed in most situations. However, a few comparisons — 
were made using the Darcy Weisbach equation. The data included 
term appropriate for the head loss expression employed. 


Computer programs were developed for each of the five algorithms studies. 
These programs were designed to solve the basic loop and node equations using — 
the | procedures previously described. Identical input data was used for all 4 
algorithms. The programs included algorithms to formulate the equations (deter- 7 
mine pipes connecting junctions, forming loops, and connecting fixed grade a 
nodes) from the connecting node data for the pipes. The programs also generated 
the necessary initial conditions. This required a balanced initial flowrate for 7 
the single and simultaneous path methods. An algorithm was developed that | 

was designed to generate an arbitrary but reasonable initial balanced flowrate 
for these two methods. For the linear method an initial flowrate is needed Jf 
but was not required to satisfy continuity. For this method an initial flowrate Jf 
_ based on a constant flow velocity in an arbitrary direction was assigned. For 
‘Single at and ‘Simultaneous methods some initial values for > 


| | 
| 
| 
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head are required to initiate the solution, and an algorithm was developed to 
provide arbitrary but reasonable initial values. It can be presumed, however, 
that calculations were often initiated with values of flowrate or head which 
are significantly different from the correct results. 
¥ Convergence Criterion. —For each method, the trials are continued until the a 


specified convergence criterion was met and the same convergence criterion ay 
5 applied for each method. Since updated flowrates are easily computed | 


for each method the change in flowrate between successive trials was used 
to check convergence. The specific convergence criterion employed is referred 


to herein as the relative accuracy and 


_ Here Q = the flowrate obtained for a trial; and Q, = the initial flowrate 
used ¥ which was computed from the previous trial. This criterion | roughly States 
that when the average change in flowrates between successive trials is less 
than 0.5% the calculations cease. This convergence criterion is more stringent 
than ones normally applied in practice. It does not, however, assure that the 
flowrates are within 0.5% of the correct values. Numerous other criterion can 
__ be applied to determine the acceptability of a solution and others have suggested 
7 and employed. Additional description of convergence criterion will be presented 
.. Accuracy of Solutions.—Each of the methods for analyzing pipe networks — 
yield approximate solutions. A solution is accurate when all the basic 1 snore 
are satisfied to a high degree of accuracy. The mass continuity equations are 
_ exactly satisfied for the three methods based on the loop equations. Each of ; 


these methods are then designed to satisfy the energy equations through an 


__ iterative procedure and the unbalanced heads for the energy equations is evidence = 


equations will be assumed to be one where the average unbalanced head for 
_ the energy equations is less than 0.01 ft (0.00328 m). For the methods based _ 
on the node equation iterations are carried out to Satisfy mass continuity . 
; junction nodes and the unbalance in continuity is the most significant indication — 
All five methods were compared to an exact solution for sy systems of under 
100 pipes. This involved a total of 60 comparisons and the results for each — 
onan were tabulated and compared as depicted in Table 1. These are 
7 _the results for the 14-pipe system which is depicted in Fig. 1. The exact solution 
was obtained in every case by carrying out the linear method one additional 
_ trial after the relative accuracy of 0.005 was reached. This table summarizes 
information which is essential in evaluating the effectiveness of each algorithm. 
The average flowrate and head range for this solution is output. Next the flowrates _ 
and heads at junction nodes are compared to exact solutions and average and 
maximum differences given for each method. The percent average and percent 


TABLE 1.—Comparisons of Flowrate and Grades for “4 pe System a m = 3. :. 


ft, 1 m*/s = 1,000 L/s = 35.3 ft/s) 


_ Average dif- 
ferences 
Percent average 
Maximum dif- wh novi 
ferences 
Percent maxi- 
mum dif- 
ferences 


Heads (in meters) 


Averages dif- 

rat 
Maximum dif- 
ferences bs 

Average 
diff / head 
 range-in per- 
— 

in 


£,.% 


number] flows | Linear | ence |SPATH| ence | PATH | ence ence ence 
— 
273.35] 273.35] 0.0 | 273.36] 0.01 | 273.31] 0.04 | 273.36] 0.01 | 272.60/ 0.75 
149.21) 149.21] 0.0 | 149.21] 0.0 | 149.25] 0.04 149.22] 0.01 149.16] 0.05 
36.21] 0.0 | 36.21] 0.0 | 36.25] 0.04 | 36.22] 0.01 | 36.72] 0.51 
| 1.66] 1.66] 0.0 | 1.66] 0.0 | 1.71] 0.05 | -4.04) 5.70 | -6.87 
-$5.34) 0.0 ~55.34] 0.0 | -55.29] 0.05 | -55.25] 0.09 | -56.40] 1.06 
95.23) 95.23] 0.0 | 95.23] 0.0 | 95.07] 0.16 | 95.18] 0.05 | 95.68] 0.45 
0.0 |—139.24) 0.03 |—139.24] 0.03 |-141.01] 1.80 
0.0 |-110.29| 0.0 }-110.26] 0.03 |-110.28] 0.01 |-110.54) 0.25 
9 =| 258.24) 258.24] 0.0 | 258.24] 0.0 258.23| 0.01 | 258.24) 0.0 | 258.44] 0.20 
124.15} 124.15] 0.0 124.15] 0.0 | 124.06] 0.09 | 124.14] 0.01 | 124.09} 0.06 
| 60.76] 60.76] 0.0 | 60.76] 0.0 | 60.79) 0.03 60.72 0.04 | 67.61 685 
12 | —17.11) —17.12] 0.01 | -17.11] 0.0 —17.22] 0.11 | —17.08] 0.03 | —15.78] 1.33 
| 531.59} 531.59] 0.0 | 531.59] 0.0 | 531.54) 0.05 | 531.59] 0.0 | 530.85) 0.74 
90.95} 90.95} 0.0 | 90.95) 0.0 90.97} 0.02 | 90.95) 0.0 | 92.30) 135 
0.00 0.05 0.43 
on | flow} flan! fos 
| 
61.49] 61.49 0.0 | 61.49 0.0 | 61.51] 0.02 | 0.0 | 61.73} 024 
: = 2 | 40.58} 40.58] 0.0 | 40.57] 0.01 | 40.61] 0.03 | 40.58] 0.0 40.84] 0.26 
31.86] 31.86 0.0 | 31.86] 0.0 | 31.66] 0.0 | 31.86] 0.0 32.06 0.20 
| 31.33] 31.33 0.0 | 31.33] 0.0 | 31.33) 0.0 | 31.33) 0.0 31.47) 0.14 
31.33) 31.33] 0.0 | 31.33] 0.0 | 31.33] 0.0 | 31.34) 0.01 | 31.45] 012 
| 36.44 36.44) 0.0 | 36.44) 0.0 | 36.44) 0.0 | 36.44) 0.0 | 36.64) 0.20 
0.0 | 32.41) 0.0 | 32.40] 0.01 | 32.42] 0.01 | 32.53) 0.12 
41.42 41.42] 0.0 | 41.42] 0.0 | 41.42) 0.0 | 41.58) 0.16 
od? 000} oo} | ow 
foo | 0.03  |00} | 086 


>, 


AG. 1 —Example iadkeie a m = 1,000 mm = 3. 28 ft, 1,000 L/s: = 35. 3 Ve) 
menionens differences are based on the average flowrate and head range and 
F these values are more useful for relative comparisons. Following this, be tne] 
data i is given in Table 2 which further relates to the accuracy of these solutions. : 
; The unbalanced heads for each of the six equations for this example are 
summarized for the four solutions (including the exact one) obtained using a 
4 ES loop equations. Maximum and average values are also given. For the two solutions ey 
7 based on the node equations the unbalanced flows at the eight junction nodes : 
are summarized along with maximum and average | values. Finally the number ty 
Ba of trials required and the accuracy attained for each of the six solutions is 
summarized. All the solutions for this example are quite good which is indicated _ 
_ by the good comparisons with the exact solution and also the small unbalanced | 
heads and unbalanced flows obtained. The worst solution was attained with "ir 
ed the single node adjustment method for which the average error in flowrate 4 


4 


only 1.23% and in grade was s only 0. 
_ Similar comparisons were made using a relative accuracy of 0. 00s for ak 
situations and the results are summarized in Table 3. Since attainment of the 
convergence criterion is not assured some liberal upper limits on the number — 
_ Of trials allowed were imposed and these are noted. Calculations were terminated . 
_ when the accuracy of 0.005 was attained or the limit « on the number of trials 
- The solutions were considered to compare favorably with the exact ape ON 
ie if the average percent deviation based on average flowrate and the head range ~ 
_ from the correct solution for flowrates and heads did not exceed 10% and me 
the maximum percent deviations did not exceed 30%. A failure was considered _ 


have occurred if the specified relative accuracy is reached and those conditions _ 
ba are not met or the maximum number of trials are run without attaining the y “ 
‘specified accuracy. number of for each are in 
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TABLE 2.—Checks on Accuracy and Convergence—14 Pipe 


Unbalanced heads | (in meters) | 


Ata 


§ 


Unbalanced fi flows (in liters per second) 


1.4 
method 
3.02 
«ape 
78 


stony: ebotom 


| 
a 
| 
= 
ind: 
Maximums 
total 


1981 


a relative. accuracy ‘of 0.005 was exact. “The number of trials 
_ not depend on the size of the system and averaged around six for the 60 - 

‘The S-PATH method elec has excellent convergence characteristics and only | 
one failure occurred. For this case a constant power pump was presented which | 
operated on a very steep head-discharge curve. . This was a low horsepower — 

- pump operating at a very low discharge. The steep ep gradient caused convergence v 
_ problems for all but the LINEAR method. This convergence problem probably — 
would not occur for the S-PATH method if the pump operated on a flatter + 
_ head-discharge curve. . Additional trials will not result in the attainment of an 

a The larger systems of greater than 100 pipes were analyzed using the LINEAR 

. ; and S-PATH methods and all attained a good solution in a reasonable number 
of trials. For the 31 situations analyzed the LINEAR method required an average — 
of 6.4 trials and 5. 4 sec and the S-PATH method required an average of 8.5 
trials and 6.9 sec. This shows that both methods obtain accurate results: with 
few ew trials and use comparable computer pte 


Accuracy = 0008 


Trials 
attained | failures | allowed| runs — 


as 


_ The three other methods studied exhibited significant convergence problems, 
and the frequency « of f problems inc increased a as + systems \ were analyzed. Since : 


A methods using systems of less than 100 pipes, results for the larger systems 
_ were not compared for the PATH, NODE, and S-NODE methods. A number 
of these larger systems were analyzed with these methods and a mgnticnnt 
number of convergence problems were encountered. However, excessive compu- 
ter times would be required to” completely “document these results so the 
. a comparisons for these methods were limited to the smaller systems. ft 
_ The eight failures obtained for the PATH method all reached the specified > 
wa In addition, for six of the pipe systems the solutions attained a relatively = 
low average unbalanced head of 0.37 ft (0.11 m) for the energy equations. 
Yet significant errors exist in these solutions. These eight cases were rerun 
with relative accuracy of 0.0005 specified, and with additional trials, all attained = 
this. The results are summarized in Table 3. Five of the cases improved to 
an acceptable degree. Two of the three failures are characterized by large _ 
unbalanced heads. However, the one additional Rae. was documented fc for 


| Trials Accuracy| Number 
‘Method allowed attained | failures 
4s 


= head \ was as only 0 .06 ft (0.02 m), yet some flowrates \ were in — 
error. This was especially evident for flow into and out of storage tanks ae 


these results are presented in Table 5. Since this prediction is of significant — 


engineering importance the results from the PATH method are reer od 4 
a This failure occurs in 1 spite of the very high aren of relative accuracy | attained =a 


- distribution } systems and yet the PATH method encounters difficulty with and 


_ When the S-NODE method is successful a a highly accurate solution is obtained > 


in relatively few trials. However, a significant number of failures were a “4 
in Table 3 for this method. The patterns of failures is the most varied of the = 


methods investigated. Six of the 18 failures noted attained an extremely accurate 
calculation of the heads. In these cases the calculation of most flowrates were 
also accurate. However, a few exhibit large errors. This | occurs when lines 
_ which carry significant flows at very low head losses are present. This situation — 
‘'s not uncommon in many water distribution systems and small errors in head 
calculations can result in a few highly inaccurate calculations of flowrate. A 
second type of failure is one where oscillations occur and the specified accuracy . 
is not attained regardless of the number of trials carried out. Table 4 summarizes _ 
results where ‘many additional trials were carried c out it for the S- NODE ‘method = 


noted for this n method, and previously published peferences to this type of " 


TABLE 4 4.—Summary of Results—Additional Trials 


‘Trials = | ~Number | Accuracy | Number 


TABLE 5.—Results for Tank Flows—79 Pipe System—Fiows in gallons per minute 7 


1 = 35.3 ft®/sec, 1 ft?/sec = 448.9 gal/min) 


Accuracy 
Number of trials | 
Flow from jom source x 813.5IN 


Tank B 897.4 


= 
Correct solution So Solution fr from PATH method 


q 84.1 


cuneunesi difficulty were cited (3,13,16). For two systems, failures of this 


type occurred for several situations at the same time that a very accurate solution 
was obtained for other situations. For one system a convergence failure occurred — 
: for the first situation. For the second situation where a single additional inflow _ 
representing a new storage tank was added and an excellent solution was obtained _ 
: in four trials which satisfied the required relative accuracy and had an average 


error in flowrate of only 0.18%. A third situation where this inflow was reduced 
a by half v was analyzed Starting with very good values of! heads $ (those just : obtained) 


similar situation occurred with another system where a very accurate solution 
was obtained for the second of three situations while the other two failed. 
In this case the third situation introduced a pump which significantly es 
the system. These results document the very sensitive behavior of the algorithm 
Failures for the S-NODE method are evidenced by relatively | high values for 
unbalanced flows at junctions, and this is probably the most reliable indicator 
- Based on the criterion for a successful solution the NODE method was as the 
least reliable of the methods studied. Although the specified accuracy was reached 
in most cases, the solution failed to meet the required standards. nd is not 
: summarized | in Table 3 but the solutions which failed often 1 reached 


which failed for three situations reached an average head change of 0.03 ft 
0. 01 m) along with a relative accuracy of less than 0.005, and yet the solutions 
contained large errors in semen In many of the situations oe as failures 


average head change for successive trials is a reliable indicator of an acceptable 
_ solution for this method. The relative average and maximum unbalanced flowrates nt 
the junction nodes are more reliable indications.§ = 
For the NODE method some selected systems which reached the specified 
_. of 0.005 in a reasonable number of trials were run at an accuracy 
of 0.0005, and for all but one system (two situations) a good solution was 
attained. These results are summarized in Table 3. It appears that most of 
the solutions which failed at a relative accuracy of 0.005 but reached that accuracy 
in a reasonable number of trials could be improved to provide an acceptable 
‘Tesuk if a a more stringent convergence criterion is met. , Convergence may be 


_ the solution even if a more stringent relative accuracy could be met. 

4 summarizes a few situations where up to 1,000 trials were carried out with 
very poor results. This implies that the NODE method is not capable of solving 
these systems regardless of the number of trials carried out. It appears that 

most significant factor adversely affecting of the NODE 


is the presence of low resistance lines. 


Significant convergence problems were documented the P PATH, ‘ale 


| 
a 
x 
i 
- = low resistance lines carrying significant flows are the principal cause "1 | 
? q 
= 
MOWEVE!, and it may Oc OVIde assurar that the 
| ’ is a good one. In the cases where the relative accuracy of 0.005 was not s 


and S-NODE methods. These methods are widely and of this 
_ Study indicate that great care must be exercised when employing these methods. | 


2 The single adjustment methods (PATH and NODE), which are usually employed 


if a hand solution or a small computer is used, must be carried out to much : 
greater accuracy than normally called for to increase e the probability that a 
methods greatly reduced the chances of failure. Stringent convergence require- 
ments for unbalanced heads for the PATH method and unbalanced flows for q 
the NODE method also may be employed to improve reliability. However, 
the attainment of a stringent convergence criterion does not assure that the 
solution is accurate and a significant number of situations were documented 4 
for both methods where substantial and serious errors occurred after ‘mecting 
Stringent convergence criterion. A few situations were documented where the — 
specified accuracies could not be attained and the solutions were not acceptable. 
De It is concluded that if a specified stringent convergence criterion cannot be | 
met using single techniques, the solution is not reliable. 


7 good solution: is attained. Attainment of a relative accuracy of 0. 0005 with these 


4 reasonable convergence - criterion n, and if this occurs in a limited number of _ 

_ trials (40 in this study) additional trials are usually of no benefit. The failure % 

rate was quite high with this method and the use of results obtained using re 

& method is not recommended unless a good accuracy is attained in a reasonable 7 
. number of trials. The best indication of an acceptable solution appears to be ¥ 

the average relative unbalanced flow at ‘the junction nodes. The oe 
are that this value should be less than2%. 

_ Each of the three methods which significant convergence 

requires a set of flowrates or heads to initiate the solution and the chance P 7 

of failure can be reduced if initial values are employed which are closer to. ty 


initial conditions does not assure convergence. i as evidenced by failures for ‘the: 
S-NODE method for situations starting with solutions for similar situations. — ; 
In some cases where line losses vary greatly or pumps operate on steep curves - 
_ the correct solution cannot be attained unless the initial values are very close 
to the correct ones and this clearly is not possible. Many of the ‘reliability > » 
problems documented using algorithms based on the node equations ee 
_ and S-NODE methods) are due to the inability of these methods to handle 
_ low resistance lines. For these lines small errors in head calculations may produce 


: ‘is inherent for r algorithms based | on node equations. This is due to the fact 


that solution 1 algorithms for these equations do not incorporate an exact | continuity 


_ It is possible that convergence can be improved and failure rates decreased 
‘ for the PATH, NODE, and S-NODE methods by introducing considerations 
such as those which have been previously suggested. In this $ study | the algorithms — 
_ for pipe network analysis were applied in a basic manner. tio 
=f given to the various techniques for improving convergence. It does appear that 
_ many of the proposed techniques are deficient. In any case the technique pel 
be thoroughly tested if the reliability isto be established. 2 —™ 
Both the S- PATH and LINEAR methods 


ut 


q | 

| 
g 

| 


lowrates and heads are computed in a few trials with great accuracy and 
the attainment of a relative flowrate accuracy of 0.005 is adequate to assure 
this. There is, however, no absolute assurance of convergence. Many piping — 
systems incorporate features which increase convergence difficulties. Since 
_ gradient methods are used to handle nonlinear terms, convergence problems 
are always a possibility, particularly if ill conditioned data such as a very wide > 
range of pipe line resistance or poor pump descriptions a are employed. However, 4 4 
this study indicates that the occurr se 
methods is unlikely. It is concluded ‘that, if possible, either the S- PATH or 
a LINEAR method should be employed for pipe network analysis and that _ 
convergence is virtually assured if reasonable datais employed. 
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= of fixed grade n nodes; 
head at system node; 4 
energy loss in pipe; 


a 
= _ number of loops in pipe system; 
- exponent for head loss expression; 
number of pipes joining junction node; 
= number 0 of pumps in pipe 
= number of pipes in pipe system; 
= in pipe; 
ow leaving system at junction node; 
= net flow into junction node; nae 


howe 
net flow leaving junction node; and by 


head difference for energy equation aiieieii's fixed grote nodes. nodes. 


a,b = junction nodes in pipe system. 
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BounDaRY (CALCULATIONS ¢ OF SLUICE 

Alexander H-D. Cheng,’ James A. Liggett,? M. ASCE, 

and Philip L-F. Liu,’ A.M.ASCE 


lize classical hydraulics problems reapp 
-more efficient or accurate (or both) solutions are discovered. Because of ‘ecie 
inherent difficulty, gravity driven (as opposed to purely pressure driven) free-sur- " 
4 — face flows have never been solved as accurately and efficiently ; as is s necessary — 
for so some design ‘and measurement "purposes. The difficulties, as pointed out 
: by von Karman (21), arise not only from the nonlinear character of the boundary 
- condition, but also from the fact that the boundary is not known a priori. scl 
_ One of the major achievements of classical hydrodynamics is the use of | 
conformal transformation in solving a simplified version of the problem. If 
x. - Bravity is absent and d the { fixed part of the boundary is formed by plane s surfaces, 
n the orthodox treatment of G. Kirchhoff, H. Helmholtz, Lord Rayleigh, 1, and 
= is able to linearize the boundary condition. However, if the fixed part 
_ of the boundary is curved or the gravity effect is not neglected, one is again 
confronted with a nonlinear problem. ~The mathematical difficulty is so far 
; insurmountable if both effects are to be fully considered. It seems clear al q 
under such conditions one must resort to numerical methods for the solution. _ 
7 Serious numerical solutions to free-surface flow began with Southwell and 
Vaisey (18), which involved a heroic effort using finite differences with hand 
a calculation. Other finite difference solutions include those by McNown, ef al. — a 
(16), and Cassidy (3). Finite element solution began with McCorquodale pore " 
‘ Li (15) for sluice gates. Many other papers appear in the literature including a 
those by Larock (12), Isaacs (11), Ikegawa and “Washizu (9), Chan, et al. (4), 
- Taylor, et al. (19), Betts (2), Diersch, et al. (6), and Varoglu and Finn (20). — a 
Ss "Grad. Student, School of Civ. and Environmental ee Cornell Univ., Hollister *) | 


?Prof., School of Civ. and Environmental Engrg., Co Cornell Univ. Hollister r Hall, Ithaca, Hthaca, 
* Assoc. Prof., School of Civ. . onl Environmental Engrg., Cornell Univ., Hollister Hall, : 
- Note.—Discussion open until March 1, 1982. To extend the closing date one month, - 
a written request must be filed with the Manager of Technical and Professional Publications, 
: ASCE. Manuscript was submitted for review for possible publication on September 26, " 7 
_ 1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of the 
American Society of Civil Engineers, © ASCE, Vol. 107, No. HY10, October, 1981. 
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‘Since flows of this nature are “characterized by highly curved boundaries, 


boundary condition is present, the Boundary Integral Equation Method (BIEM) 
has proved to be a superior numerical technique (13). Direct application of 
_ the current methodology to transient water wave problems poses greater numerical 7 


av  Miotky' in terms of instability due to the additional nonlinearity from a dynamic 
hyd surface condition and the fact that there is no natural damping. For many 
hydraulic problems in open-channel flow, particularly those associated with flow 


3 measuring devices, one is concerned only with the steady-state properties, such 


as the steady-state free surface profile and the head-discharge relationship. In 

‘such cases it seems unnecessary to solve the transient problem. Instead, a 
Scheme based on a perturbation principle | is employed in this investigation to 

provide rapid convergence of the solution to problems of sluice gates and 2 
_ spillways. The BIEM, at least as now applied, is not the final answer to these 
classical problems. It does, however, represent a further step toward their solution © 
the of physical models with numerical calculations. ‘sa 

du ti 

In the sluice gate and spillway problems only the ‘steady state i is considered. | 


which = ‘Laplecian operator; and = the stream function. Both 


- the fixed boundaries and the free surface are streamlines, , therefore w is taken 


on the free surface and on the sluice ( 


in which Q = the flow rate } per wi unit wi width. On the free > surface the dynamic. 


which v = velocity; g = the of f gravity; 


=, 


in which n = the unit normal from the free surface; > Eq. 3 becomes" 

SSS 


A 
either a fine mesh or the use of higher order interpolation functions is usua 
required for a numerical approximation. In addition, due to the intrinsic nonlin- _ 
aq earity of the problems, a formidable number of iterations seemstobe unavoidable 7m 
_ in the solution procedure. Thus, the efficiency of the solution method is important : 
: and can become a major factor in the selection of anumerical method. ==  ####$#7m 
For many problems in porous media flow, especially those involved in the oa - 
= determination of unknown location of a free surface where a nonlinear kinematic © o — 
| 
— 
| 
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BOUNDARY CALCULATIONS 
the suriace 
Either the flow rate, Q, or the Bernoulli constant, B, are known; the other 
Both the "sluice gate and ‘spillway flows. extend to +o. For the purposes 
of the numerical solution the inflow and outflow streams are cut at right angles 
to the primary velocity. On the cut the condition 
: © 


is  eatiel which ‘means that there is no velocity normal oe main flow. This — 
7 condition, although approximate, is applied sufficiently far from the sluice gate 
: or the spillway crest that it is highly accurate and any small error does not 
affect the interesting part of the flow. 
A description of the use of the BIEM using Eq. | with boundary conditions — 
of the Dirichlet type (Eq. 2) and the Neumann type (Eq. 6) is found in Liggett — + 
(13). The addition of the free-surface boundary condition (Eq. 5) forms an — 
: "additional complication. This factor is treated in much the same way as it has — 
been for most of the finite difference and finite element solutions. That is, 
for problems with a known flow rate, Q, the position of the free-surface boundary _ 
_is assumed and the problem solved from Eqs. 1, 2, and 6. The “constant,” — 
B, of Eq. 3 is then calculated on the free surface. if B is the same for all ~ 
_ free-surface points, the f problem i is solved. Otherwise, ‘the assumed free. surface ae 
_ is adjusted, iteratively, in order that B becomes constant at all points. The 7 
_ method of adjusting the free surface is problem dependent: for the sluice gate 
the adjustment technique can be made efficient; on the spillway a 1 


is part of the solstion, If Qisa constant for all points, the 
"problem is solved; otherwise, an iteration scheme must be used to adjust the 


Solution Procedure.—The geometry of a vertical sluice gate is shown i in Fig. 
1. If the flow rate is known, the iterative scheme for adjusting the free- surface 


 Jocation can be described as follows. Introducing v as the average 


at a vertical cross section, the flow rate can be expressed as 


For the purpose of the free surface adjustment, iti is assumed that 9 is ee 
free-surface velocity, v. (Although this assumption i is poor in the neighborhood 
of the stagnation point, S, of Fig. 1, the pertinent term in the iteration equation 
: is so small in that neighborhood that it does not spoil the convergence.) A 


Eq. 3¢ can be w written in terms of the depth 
yo 


| 
_ constant, B. The problem is solved first by assuming the location of free surface § | 
: 
| 
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4 in which terms i in ‘eqnr and higher order have been neglected. lected. Subtracting 


4 Eq. 8 from from Eq. 9 and solving for Ay 


hag: 


oo 


FIG. —Definition Sketch of Vertical Sluice 
of the gate and greater than unity on the downstream side. Sad 


P 

buted 


We 


1. The position of the free surface is estimated and ae are "distri 
on the boundary, including the free surface, papiirstad 
2 The BIEM is used to solve for v = dW/dn at each rane surface point. "ee 3 
3. Using Eq. 3, the Bernoulli constant, B,, i is computed for r every free-surface * 
point, i, including the tip of the sluice gate, point C (Fig. 1). Since point C 
cannot be adjusted, i. €., y = bis fixed, yet Eq. 3 applies at C, it is assumed, = 
temporarily, that B* ' (where the superscript indicates - iteration number) is 


the correct value of. Beverywhere forthe flow. 
-. Using Eq. 10 the i-th node on the free surface is displaced an amount 


When |AB*“ which € is some small number, the iteration process 


considered to > hav v ged; otherwise, 24 are repented. 


| 


Po assumes that moving one free-surface point does not freon the Berouli 

constant at other points. Such an assumption is incorrect since all points are 

_ A similar method is used to solve for the flow rate given a Bernoulli constant. _ 

d The flow problem is solved by replacing the kinematic free- surface « condition, ; 
Eq. 2, with the dynamic free-surface condition, Eq. 5. The solution yields condition, 
values of the stream function at each of the free-surface nodes. A temporarily — 
“‘correct”’ flow rate, Q**', is needed for performing the free-surface adjustment. 
_ As in the previous calculation point C is used and in this case the stream 

function at point C, y*, can be a good estimation of the correct \ value; however, 4 

a better rate of convergence is achieved by taking the wie ; 


factor which can be used to accelerate the rate of convergence, must be determined 
i empirically. (Eq. ul can also be written — in terms of the Bernoulli constant 
‘to improve the prediction of B“*' for the case that the flow rate is given and 


The free-surface adjustment scheme, which to Eq. for the 


st 


= 
jawiich | Ay A 
Linearizing and solving for Ay, produces = 
107) 
For large gate aguaiuas,. hawt, the Froude number for portions of the 
_ downstream profile is not much larger than unity. In those cases the Ay, of 
Eq 14 tends to be too large. A more conservative , yet si stable, scheme i is used 


Vv re: 

. ‘Thus, E Eq. 14 or Eq. 15_ leads to a new location of the free surface in the 

same way as Eq. 10. The remainder of the procedure is ; similar except that 

accuracy checks are done on the flow rate instead of the Bernoulli constant. 
Though the aforementioned methods give a reasonable rate of convergence, t 

it is found that the stream function and the free-surface velocity (when obtained 


as a result of the calculation) are sensitive to the free-surface location in the _ 


subcritical and the supercritical flow re regions, , respectively. The assumption of | 
a small perturbation, as used to derive Eqs. 0, 14, and 15, may fail if the 


4 
| 
&g 
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= initial guess of the free- surface location is too far from the true solution. The 
_ process is improved by applying each method only to the region where it is 
most efficient. This method, which calls for a mixed boundary condition « on 
the free surface, is presented in the following. GESTS TESA 
a Eq. 2 is applied to the upstream free surface to give a Dirichlet type boundary — 
q condition and Eq. 5 to the downstream for a Neumann condition. Since — 
; one of the two quantities, Q and B, as required in both Eq. 2 and Eq. 5, 


- is known for a given problem, the other is replaced by a temporarily ‘‘correct” _ 
value, , determined from the aforementioned equations. The free- 


initial profile is a 
= quarter of an 


» Stag 


«AG. 2. —Convergence of Downstream Profile for Vertical Sluice Gate 
surface adjustment schemes, Eq. 10 upstream of the gate and Eq. 14 or 15 


downstream of the gate, are used to find the new free-surface location. Both 
SS of sluice gate problems, those with known flow rate and those with st 
b Bernoulli constant, can be treated by this method. A rapid rate of convergence 
< is always observed, even with an initial guess which differs greatly from the © ud 
- All the previously mentioned solution procedures require an initial guess of 


either the f flow rate or the Bernoulli constant. Although equations similar to ad 


Eq. 11 can be used to give a temporarily ‘ “correct” value for each of the 


g 
— 
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BOUNDARY CALCULATIONS 


subsequent i iterations, an alternate procedure is necessary for the initial iteration. 


An empirical equation will suffice but another starting technique has been used © 
in the calculations presented herein. For the case of known Bernoulli constant q 
the boundary conditions for the first iteration are modified as follows: af/an 
computed from Eq. 5 is ‘specified on both free surfaces, AS and CD. On _ 
sluice gate, SC, d/dn is prescribed using an arbitrary interpolation formula, 
e.g., linear or quadratic, between the two known values at S and C. Thus 
a well-posed problem is defined which does not require an arbitrary guess of “ : 
an initial trial value of Q. oa problem with given flow rate can be treated 


= Sluice Gate Solution. —To compare e the accuracy and efficiency of the BIEM 
"solution with a most recent FEM solution (11), the following problem is solved. 
TABLE 1. —Downstream Profile of Vertical Sluice Gate with b/B = 0. 3 


in 
_ Squere feet 
-~persecond 


0.278 1.304 


0.031 193 1.303 


1.303 


1.303 


0.179 


Isaacs 
"Solution of the 6th iteration. 
Note: as 1 ft; Q (Isaacs) = 


The atee gate geometry is given as shown in Fig. 1 with B =1 (arbitrary 
length unit) and b =| 0.3. The initial trial upstream free-surface p profile is flat i 
(y = 1) and the downstream profile is fit by a quarter of an ellipse with an 

arbitrary downstream depth of y = 0.56 (Fig. 2). The optimum value of a, 
as defined in Eq. 11, for each problem can be found on a trial-and-error basis. 
For all cases of sluice gate flows solved in this paper, a = 0 is found generally 
Satisfactory, though v values ranging from 0-2 tend to speed the rate of convergence. i . 
_ ‘The flow boundary is discretized into 39 linear isoparametric : elements. The 
problem has also been solved using cubic spline elements as described by Liggett — 
and Salmon (14). Since the results using the different elements do not differ, 


only those using linear elements are presented herein. gr 190, en . 


yin in square feet 
0.000000 
| om 
0.139618 0.209 | 0.208 07 
0.300813 (0.187 7 0.186 1,305 
| 0180 | 0.180 1.304 133 
20000 | | 1302 
q 1.304 cfs/ft; and Q (BIEM) = 1.303 cfs/ft; Qis the {Jf 
— 


rather arbitrary profile, is achieved at the 3rd iteration 

qi e., free surface is adjusted only twice). Though not shown, the upstream 

profile converges to the solution even more quickly. 
>> computation is carried to the 6th iteration where the error, € = max _ 


yi / Q**"), is within +0.1%. The results for the downstream profile are given 


eel FIG. 3. —Upstream Free- Surtace Profile for Sluice Gate a 


oh 


aah 


FIG. 4.—(a) Pressure Distribution along Channel ae for § Sluice Gate; ; (b) Pressure 


Distribution slong Stu Sluice ot od aks ony at ai Thay 


* in Table 1. Also listed in Table 1 are results from Isaacs’ (11) FEM solution. 


The FEM solution is obtained using approx 150 elements, some of which are 
_ curved cubic triangular elements (10), and 30-40 iterations. A similar finite element 
= by Diersch, er al. (6) used 120-125 quadrilateral elements, corresponding 

to a maximum of 1,061 degrees-of-freedom, and 17-25 iterations. As demonstrated a 


™ the efficiency of the BIEM solution i in terms of the mamber of eee, 
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| | 
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BOUNDARY CALCULATIONS 
number of elements, and computation time is ‘clearly superior to > the FEM o- 
comparable accuracy. The outstanding accuracy for | such low order elements — 
_ may be attributed to the fact that all numerical approximation is confined to 
the boundary and the integration is exact. The upstream free- surface profile 


Profile of Radial Sluice Gate 


fot |  BIEM 
| 


154893] 2.51 


2.04077 g 
3. 51683. 2. 7 
2.09 
$.03810 2. 6 
5.42010 | 2.06 


Solution of 8th iteration. 
a Note:  s = Bos 556 ft; Q (Isaacs) = 


Fongmeier a Strelkoff 
BIEM 


b/B=6 


ynstream Free-Surface Profiles for Various Gate Opening Ratios s 
C. 


and the geessure distribution along the sluice gate and the bottom are plotted 
in Figs. 3 and 4, respectively. va with 


A radial s sluice gate fi flow problem, Solved by Isaacs (11), is also Presented 
in Isaacs’ " Paper. Thirty-nine li 


| 
@ 
41.96 cfs/ft; and Q (BIEM) = 41.92 cfs/ft; Q- 
4 
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profile of the 8th iteration (« = +0.05%) is compared w 
FEM solution in Table 2. The use of linear elements to , adi a highly 
curved geometry isindeed adequate. 
Vertical sluice gate flows with various gate opening ratios can also be solved. 
The downstream profiles for those with b/B ratios, 0.1, 0.2, 0.3, 0.4, 0. 53, 
and 0.6 are sketched in Fig. 5. Also shown for comparison are analytical solutions _ 
by Fangmeier ar and Strelkoff (7). The agpocment between the BIEM and analytical 
The coefficient for a ra sluice gate is 13 
in which d,, as Pom in Fig. 1, = the paar oe flow depth. The 
contraction coefficient as a function of the gate opening ratio is plotted in 
my 6 along with the analytical solution by Fangmeier and Strelkoff (7), numerical 
solutions by Southwell and Vaisey (18), and Isaacs (11), and —— or 
Southwe 18 Vai 


segments were to cathe! flow b Th 
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A FIG. 6.—Sluice Gate Contraction Coefficient Versus Various Gate Opening seed 
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Overall. themipdiway ‘solutions ave considerably ors diffiow: 
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FIG. 8.—Bevelled Tip which is used to Represent Sluice Gate with Partially ue 
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FIG. 7.—Sluice Gate Discharge Coefficient Versus Various Gate Of 
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by Benjamin (1) and Rajarstnam (17). The are larger 


the theoretical ones. A possible is in the following 

discharge coefficient for the sluice gate is” 


in which d, = ti upstream uniform flow depth, as ibe in Fig. 1 . Compariso 


between discharge coefficients obtained theoretically (6, 7, and 11), and = 
mentally (8) is shownin Fig.7. = 
Numerical Experiment. —Studies (il, ,12) how that the thickness 
caused by the boundary layer growth on the channel walls is not large enough — 
to account for the discrepancy in the contraction coefficient as shown in Fig. 
6. A numerical experiment is is designed to investigate the the effect | of a change 
the local geometry of the gate tip. 
A gate with a bevelled tip is sketched in Fig. 8. The distorted part is made 
ns x small so that it does not affect the silane ratio. No special treatment of singularity 
} 


- TABLE 3.—Contraction and Discharge Coefficients for Be Bovelied Gates” 


0.550 


is given to the corner where the gate is this 

a rounded corner. The flow is attached to the bevelled —_— and thus the angle | 
_ of flow release is 

The cases of b/B = 0.3, a/B = 0.003, and B = 0°, 5°, 14. >. wv, &,. 

65°, and 50° are solved. The contraction coefficients, Cos and discharge 
_ The cases with gate geometry the same as s the aforementioned while the | gate : 

“opening ratio is increased to 0.4 are also calculated. The range of variability 


for b/B = 0.3-0.4, and B = 0°-90° is shown as the shaded area in 


~The foregoing experiment is by no means a rigorous simulation of the physics 


near the gate tip. However, it does demonstrate, ‘qualitatively, the effec 
‘local geometry in the form of a a tip. banda? 


—Spuway lees 


the sluice the problem viewed as solving for the Bernoulli c constant 


| 4 

Zz 


the flow rate or vice- versa. the problem has done both 
ways, the majority of the solutions are of the former type. 
s Be the spillway solutions are considerably more difficult than t those o of 
4 “the sluice gate. As indicated in Fig. 9, there is a zone of uncertainty where 
a change in the water surface level has a negligible effect on the Bernoulli 
ms constant. Thus, a straightforward iteration procedure, as is done for the dice es 
Bate, , usually leads to a “‘solution’”’ in which the result of one iteration differs 
negligibly from the last but | where the water surface in and near the zone of | 
- uncertainty may bear little resemblence to the experimental results. The final 4 ; 

- calculated free-surface position using this technique is dependent on the initial 

a guess and on the way the iterations are carried forward. Thus, a different procedure 


The assumption that points can be adjusted independently w was abando ae 


A Newton- method used which is stated | 


Ay 


‘Zone 
uncertainty 


lows 


furs 


9.—Definition Sketch of Spillway TOS, 

(48) = = the departure of the ‘Bernoulli constant for each free- surface node 
7 from its target value. The target value of the Bernoulli constant for iteration, — 
ok + 1, is taken to be the average value of iteration, k. The matrix [dB/dy 5 
is computed by displacing, in turn, each free-surface node a small distance 
a from its assumed value (or the results of the previous iteration) and computing 
the change in the Bernoulli constant at all free-surface nodes due to the 
displacement. at that particular node. The calculation requires a complete solution 
for each free-surface node. Thus it is fortunate that the BIEM can do such © 
- calculation efficiently. Since the perturbation matrix changes little from one 
iteration to the next, it is not necessary to recalculate it each time. When it 
does have to be recomputed, it is necessary only to perturb and recompute 
= points i in | and near the zone 1¢ of ‘uncertainty. yn 


. in which {Ay} = the vector of the adjustment of all free-surface nodes; and 
| 
\ 
4 
| 


(tainty. Using cubic ‘spline elements leads to good ‘solutions in ‘the 
error in Bernoulli constant con be made as small as desired at the nodes. However, — 

¢ the error between the nodes, as calculated from the interpolated velocity and = 

position of the surface cannot be made arbitrarily small, although it is satisfactory. : 

: _ The best procedure thus far discovered is to place a number of ‘ ‘pseudo-nodes” =— 
between the free-surface nodes. Either linear or cubic spline elements are used, 7 
ben the elevation and velocity at the pseudo-nodes interpolated accordingly. 

7 - Only the real nodes are perturbed yet [dB/dy] is computed for all nodes. 


‘Ifthere are N real surface nodes and M pseudo-nodes then dy} becomes 


BIEM nodes for 40 cfs 


— Experimental curve Qs 59 cts 


distribution for Q= 
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10.—Surface Profiles for Spillway with Q = 21 cfs 40 cfs and Pressure 
Distribution on on on Spillway for Q = 40 cfs at r 
“both real and pseudo. The surface is adjusted only at the real ates As a 
- practical matter some limits or damping, or both, is applied to the {Ay} in 
_ order that the adjustments to the free surface are always reasonable. Of course, 
the departures of the Bernoulli constants cannot be made arbitrarily small since 
Eqs. 18 are never exactly satisfied. The solution is considered to have converged 
when the change in the surface position is sufficiently small. = = = 
_ Spillway Results. —Fig. 10 shows a standard WES spillway designed for 59.21 
cfs/ft (5.5 m°®/s/m) of width. The height of the spillway is h = 24.78 ft (7.55 
m). (Standard spillway designs are explained in Ref. 5.) The solid line is the 
experimental water surface. The BIEM results conform closely to the cnemene 
— . The computed discharge coefficient is C = = Q/B” = 4.05 whereas © 


ex 


7 
| 
| 
q 
q 
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BOUNDARY CALCULATION 


gives 4.03. It is difficult to make the. Bernoulli constant its target 
value on the last nodes of the spillway face since the velocity is high and 
a small velocity change makes a large difference. Of course Eq. 19 could be 
weighted in order to improve the accuracy in any region, but ‘the ; result would — 
be to shift the i inaccuracy to some other area. 
7 ™ Fig. 10 also shows the surface profile for a flow of Q= 40 cfs/ft (3. 2 
m’*/s/m) over the same spillway. The computed discharge coefficient is 4.05 _ 
in this case. The figure shows the pressure on the spillway which is everywhere 
‘4 positive but at a low value for the node at x = 3.5 ft (1.07 m). The free- surface 
y pressure v was also computed. Although it should be identically zero (atmospheric), 
_ it had an average value of 2.6 psf (0.0012 atmosphere) with a maximum of — 
16 psf (0. ‘Thus, Eq. 3 has been to a 
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The following symbols are used in this paper: ’ 
= discharge coefficient for or spillway; work 
contraction coefficient for sluice 
Ca = discharge coefficient for sluice gate; 
= upstream uniform flow depth; 
d, = downstream uniform flow depth; vig . 


= 


unit outward | norma ay 


Orv 


pressure; Thies 


flow rate per unit width; 


angle of bevelled sluice gate tip; and 
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INJECTIONS FOR OptimuM 


By D. Fitzgerald," A. M. ASCE and Edward R. ‘Mz ASCE 


The mixing of a miscible substance water flowing i ina is often 
part of a chemical, or mr biological water treatment t process, or r both, and i: is also 


mixing g devices such as batch r mixers, s, fixed blades in a pipe line, pumps, _ orifices, iy, 
valves, etc. are used in many operations. However, their use may require a 
- Gareption of the flow or cause a substantial head loss. Thus, for some ae. 
investigation of alternative mixing methods may be warranted. One viable 
alternative is to use the pipeline as a mixing chamber by injecting the ‘substance — 
into the pipe f flow. The method is sometimes referred to as pipe- flow mixing, ’ 
_ In situations in which enough pipe neath, is available and the speed with a 
~ which mixing is eccompliched is not important, the ambient turbulence in a 
pipe flow can ultimately provide complete mixing of miscible substances. 
Information is available on the distance ‘Tequired to achieve a certain degree 
“of mixing in ‘such situations (6,9). This { paper is concerned with those situations 
_in which more rapid mixing is required and, more specifically, with the use - 
_ of jet injections to provide initial mixing, thereby hastening the mixing process. n 
- maintained, and operated with minimal interruption of the flow. . There are no- 
concerns about structural problems or flow-induced vibrations of the injection n 


With the injection tube at the wall, a jet-injection system can be installed, 


tubes since there are no struts or tubes extending into the pipe. Also, the loss 7 


of total energy in the ambient flow due to the jet injection is usually negligible. i” 
’ In this paper, the injected fluid or additive is called a tracer. The mixing 
_ distance is defined as the flow distance from the point of injection to the point — 
: where the variation in concentration over the cross section is some : small, , Specified a 


mixing for a given | of jet injection. “Source” or “source- De-type”” 
"Sr. Grad. Engr., Turner Collie & Braden Houston, Tex. 
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refers to situations in which the injected fluid has no slapdilinne initial 1 a 

‘The objectives for this study were to: (1) Experimentally investigate several 


jet injection systems to determine their mixing characteristics; (2) identify ool 


Definitions. —A which can be used to characterize the unifor- 


mity of a concentration distri distribution over a cross section is the coefficient of — 


C= 


s 


wo 


_ in which c = tracer concentration at a point in thec cross section: t= een resctentl q 
average concentration; and A, = area of the pipe. _ The variable C, is largest 
near the point of injection | ‘where the ‘concentration distribution is highly 
‘nonuniform and decreases with increasing flow distance as the concentration 
distribution becomes progressively more uniform. For large flow distances, ; 
empirical C, values are nonzero and are a function of the random measurement 


_ errors associated with the experimental apparatus. For a given injection system 

: and given flow, the mixing distance increases as the desired C, decreases, i. = 
as the requirement for the degree of uniformity. increases. _ Both analytical and 
_ experimental evidence indicates that, for ae C,, the maximum concentration 
variation within a cross section, i.e., |c — iS linearly related toC, (9). 
_ From dimensional analysis, the 4 ne penetration, Y, of a nonbuoyant jet 

into the ambient crossflow canbe writtenas— 


meter ratio = D,/D, = . initial Mii ‘number; « a 

= yh ya of jet relative to the rie Y wall: and subscripts j and p refer 
to the jet and the pipe, respectively. In Eq. 2 it is assumed that the turbulence 4 

_ characteristics of the pipe flow do not influence the jet penetration. Experimental — 

results (1,7, 10) indicate that the behavior of jets injected into pipes is not influenced 

_ by the volume flux ratio so that Q, can be eliminated from Eq. 2. [Wright 

mony) considered the significance of the volume flux for jets in unbounded 

crossflows. ] Further, if R, is large enough so that kinematic similarity of the 

jets exists, then R, is not a significant parameter and Eq. 2 can be written 


| 
operating conditions; and (3) investigate the effects which secondary currents _ 
&g __ This paper summarizes some of the background and experimental results. 
| ant 
] | 
| | 


INJECTIONS 
Similarly, the initial mixing provided by jets is a ft function of M, and a, plus i 


‘the number and location of the jets around the pipe wall. Sources are the & 
case in which M. = 0 so ‘there i is zero penetration and zero initial 


_ by the second writer and Ger (8) with a summary given in Refs. 6 and 9. ; 
It was shown that for given locations o of passive sources, OC, is: a function - 
of a dimensionless longitudinal distance 


in which = radial diffusion coefficient, R, 


can be Ea. Sc can be as as (8) 


4 
_ in which f = friction factor; and Se, = turbulent Schmidt number. The differential — 
mass balance equation was solved analytically to investigate the diffusion of - 
a conservative, neutrally buoyant substance released continuously from point — 
sources in a steady, uniform, turbulent pipe flow (8). Except for source 
configurations which are symmetrical about the pipe center line, the asymptotic 
Slope of log C, versus Z is a function of a . 
in which e, e, = circumferential diffusion coefficient. ecieenil results were 
used to evaluate y and e,, but the values depended on the value assumed — 
for Sc,. Using the concentration data collected by Filmer and Yevdjevich (5), 
“j Clayton et al. (3), and Ger and the second writer (7), 1 was found to be 1. 8, 7 
; assuming a turbulent Schmidt number of 1.0, or a= = 1.2, assuming Sc, = 0.77. _ 
- Fig. | is a plot of log C, versus Z with y = 1.8 and Sc, , = 1.0 fora number a 
of different point source locations represented by p’ = which r = 
- radial coordinate measured from the dey center line. Thus, p’ = | represents 
_ a source at the pipe wall, and p’ = 0 is at the center line. The asymptotic 
2) slope for the =0 curve also ‘Tepresents the € asymptotic rate of 


mixing or slope of log ie versus Z for any initial concentration distribution — 
or any number of point sources provided that the initial distributions or sources - 


- are symmetrically arranged about the pipe center line. Similarly, the p’ = 1 
_ curve for large Z represents the asymptotic slope for all asymmetrical cases ~ 
and for the slowest rate of ambient mixing. For brevity, these characteristic — 
slopes will be called the ‘symmetrical’ slope and the ‘ “asymmetrical” slope. 


#s or asymmetrical initial concentration distributions or point source distributions, F 


the log c, versus Z curve can have somewhat different shapes for small Z. 


| 
Ambient Vivino —Ambient mix! 
i 
radius; and z = flow 
| 
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he curve can follow ‘the center line source eventually it ‘will break 

away to the ‘“‘asymmetrical’’ slope. See, e.g., the curve for p’ = 0.2 in Fig. 

1. The point at which the curve changes from the ‘ “symmetrical” to the > 
‘asymmetrical’’ slope depends on bent initial degree of asymmetry. — 


The * “symmetrical” Slope for p’ = 0is associated with perfect ‘symmetry, 
‘but perfect symmetry can never be obtained in experiments. Thus, it might A. 


be expected that empirical curves for log C, versus Z should eventually have . 

q the asymptotic slope for all asymmetrical conditions, but this is not necessarily . A rN 
the case for two reasons. First, Fig. 1 shows that the value of Z at which 
die curves for asymmetrical injections break away from the curve for a 


0 depends | on the value of more generally, on how nearly ‘symmetrical 


te or — condition is. For | = 0.05 and Z<0.4 (or > 0. 01), 


wed 


3 


| 


Dimensionless Longitudinal Distance, 


anne to p’ = 0.05 and pe’ = 0. Thus, a slightly asymmetrical — 
could give results which followed the ‘ “‘symmetrical’’ slope. The second —— 5 
is that at some value of C,, the empirical results begin to represent random " 
_ measurement errors” rather than degree of uniformity of the concentration — i. 


: distribution. Thus, it is normally not possible to measure extremely small values 


of C, to determine if the log C, versus Z curve has broken away from the ef 

‘ _ “symmetrical’’ slope. Empirical curves for nearly symmetrical situations can 
therefore follow the ‘‘symmetrical’’ slope throughout the range of ee > 
values. values of C, below which errors 


~ 


' or fluorescent tracers (9), 0.02 for conductivity (6), and 0. 003 for 
_ Experimental and analytical studies dealing with ambient mixing have been _ 
summarized elsewhere (6,89). 
Jet-Injection Experiments.—A turbulent jet injected with significant momentum 
the ambient flow can create rapid initial mixing and can significantly — 
4 reduce the mixing distance. Chilton and Genereaux (1) published one of the = 
- first reports addressing the use of different configurations of injection systems _ 
_ to promote rapid mixing. They used air flowing through a glass pipe and an_ 
air jet containing smoke as the tracer so visual evaluation could be made concerning 
relative mixing for each injection recognized from their 
existed 


"a Some information on the relationship of jet penetration to M, can be obtained _ 4 a 


” from Wright’s (11) work on jets in cross flows. His momentum length scale, 


Fig. 4 for the the distance from the pipe wall, indicates that 


— It is difficult oi use +e 9 to obtain an absolute evaluation of the jet penetration _ 


because of uncertainties concerning the point of transition from initial mixing © > 
to ambient mixing. Nevertheless, it is possible to use Eq. 9 to obtain relative — 7 


penetrations for different M, values by assuming, e.g., that the transition can 


values represented by the Subscripts 1 and 2. Differentiation of Eq. 9 and 


so that the relative penetration of jets into the pipe should be proportional — 
a Ger and the second writer (7) experimentally found the optimum momentum — 
‘atio, M; *, to be ) APPrOX 0. 016 for a a single c cross- “flow jet. ‘This value i is considered — 


: about me center line of the | pipe after the jet was bent over by the ambient 


, be represented by the same slope (dy/dz) of the trajectory for any two M, 


| 
_ characterizing the penetration for most of the jet injections, the trajectory can 
Gg 
| 
a _ diam galvanized steel pipe used for this study. They showed that for M —-_ ; 


flow. For M, 
the jet overshot the ‘center line. For M, _ either greater than or less than M* 


the mixing distance | increased compared to the situation in which M, a. 
_ Therefore, for a single jet-injection system to have the maximum effectiveness 
as a mixing device, it must promote rapid initial mixing and distribute the tracer 
_ symmetrically about the center line of the pipe to allow the diffusion processes” 


ene ambient mixing part ‘of ‘the process e even n with initial jet mixing 
(8). Normally, the jet is dissipated rather rapidly so that the initial mixing region 
_is rather short (only a few pipe diameters, at most, fee small diameter engi > 


ft pipe Reynolds santber was between 30,000 ond 40,000 and the friction — 


t vial 


‘TABLE 1.—Location of Sampling Stations" 4 


7. 0.053 
237. 00 ad 


approx 0.025. A nonbuoyant, Sodium chloride tracer was s injected at 


tubes with the ‘outlet. flush the pipe ‘The Gatribution of tracer 

concentration in the ambient flow was determined by conductivity measurements 4 

on two perpendicular diameters with a total of 13 points at sampling — 

spaced 10 pipe diameters apart, as given in Table 1. , NS oe 

_ The two types of injection systems investigated were: (1) A single jet originating 

at the pipe wall with the angle of injection, a, varied from 90° (cross flow) — 
- relative to the ambient flow direction to 150° (almost counter flow); and (2) 
two jets which originated at the pipe wall, were at opposite ends of the vertical 
= diameter, and had an ne of injection of 90°. After finding the optimum conditions 


motor propeller upstream of the injection point. 
a Wall Source Experiments.—Previous experiments with a single wall source 
‘ were repeated to check the experimental techniques. The data from this study _ 
_ were in good agreement with the data from the other studies (3,7) and with — 
the analytical solution for a wall source (see Fig. 1 with °’ ‘= 1.0). It was 


q 
+ 
4 
Experiments were c 
| 
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_ sampling points were sufficient to provide a an accurate saaiession of C, at 4 

of the measurement cross sections. 
Single Jet Injections.—The single jet-injection experiments consisted of an 
_ injection tube positioned at the top of the pipe with the outlet flush with the j 
if pipe wall and at angles of 90°, 120°, 135°, and 150° relative to the pipe a 


_ The fluctuations in the background salt concentration and the random errors 

oa in the measurements had a negligible effect on C, for the highly nonuniform 

; concentration distributions near the point of injection. However, the experimental 
_ inaccuracies became progressively more significant with increasing uniformity 
_ of the concentration distributions and with increasing flow distance. ed comparing 


135° Jet 
fun Symbol! H 


= 9030 
OI 
Be 


~Centerline Source 
Slope) 


Longitudinal 
2 .—Coetficient o of Variation for Single 135° Jet 
the empirical log C, versus Z curves to the general analytical behavior of log 
_C, versus Z curves (see Fig. 1), it was estimated that C, values equal to or 
less than approx 0.02 represented random measurement errors rather than degree ; 
_ of uniformity of the concentration distribution gale values of C, below a 


_ For each a, M, was varied to provide concentration data for estimating the — 

~ optimum M*. The method for selecting M* is considered at ‘the end of this — 
section. Fig. 2 shows the empirical C values obtained for a = 135°; similar 4 
variations of C, with Z were found for each a. The individual C, values have - 


been published elsehwere Fig. 3 shows the comparison of C, values for 


7 
| 
1 
4 
| a 
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as the flow distance ‘required. for C, to reach some ‘specified value. For 
equal to 0.02, Fig. 3 shows essentially no reduction in Z,, for M, = M* as 
a was increased from 90°-135° and a reduction of only about 15% ata = 
4 150°. The reason for this behavior is that transverse symmetry of the jet became 
_ progressively more difficult to maintain as a was increased. ‘The asymmetry 
could be seen directly in the concentration distributions (6), and the results” 
of the asymmetry can be seen in that the C, _ values for the larger a values F 
tend to curve away from a line with the ‘ ‘symmetrical” slope. For C, = 0.1, 
there is a consistent decrease in Z,, with increasing a since this value of C, 


“symmetrical” slope (see Fig. 3). 


o 

te 
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—Coefficient of Variation for Optimum Momentum Ratios 


_ The variable Z m for a = 150° is more than 50% smaller than for a= 


instead of a wall source. ‘Volt 
Fig. 4 depicts the variation of Z , with M, for each a and for two values - 


of C,. For C, = 0.10, Z,, is not very sensitive to small variations of M, from 
Mt, as indicated by the relatively flat curves through the data points. The : 

- indicate a greater sensitivity of Z,, to M, for C, = 0.02. The calculated — 
value of Z,, for M, = 0, i.e., a wall source, is shown in Fig. 4 for each a 
to demonstrate the reduction | = Z,, Which is _— by using jet injections. 


| 
| 
| 
| 
ect 
o 
owever, the opumum momentum rauio for a@ = is approximately six times 
: kz than that for 90°. The variable Z. for C. = 0.02 is 66% less when _ 
q 
| 


The dashed | — in Fig. 4 are very approximate, having been estimated pam : 


the information in Fig. 1, | assuming that the jet penetration is proportional t to | 
the square root of M, (Eq. 10) and that each jet can be represented by a a virtual 
source at the p’ corresponding to jet penetration. 
e _ The experimental program was designed and executed to conduct tests for 
- range of M, values which would bracket M* for each a. Unfortunately, as — 
: ‘Figs. 2 and 4 indicate, not all of the sets of experiments actually provided — 
bracketing. This is an unfortunate condition which ‘resulted from a a re- analysis 


bem 


“FIG. 4.—Mixing Distances for C, = 0.02 and 0.10 moped 
F to conduct additional tests. In some cases (primarily for a = 150°, Fig. 4) 
a clear definition of M* could be obtained from the available log C, versus 
. Z graphs. However, the values of M* for the other injection angles were selected © 
both by studying ‘the C , graphs and by checking the concentration yn distributions: 
4 at the downstream stations for symmetry with respect to the horizontal diameter 
_ since it was this symmetry which was indicative of optimum penetration of 
the jets. Due to the horseshoe or kidney-shaped concentration contours which 
develop for jets in cross flows (4), it was not uncommon for optimum conditions — 
to have the maximum concentration near the center line at Station | and below _ 


4 ‘the center line (e. g., at r/R, = 0.1) for Stations 2 and 3. ms Station 4, send 


| 
| 
| 
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oncentration were about the diameter 
a optimum conditions. Symmetry about the vertical center line related to 
pom of the jets. By thus studying the concentration distributions, all of 
_ which have been published elsewhere (6), it can be seen that M, = 0.0124 


for a a= = (9 g gave underpenetration while M, = 0.0138 gave viahinineati _ 


so that the optimum conditions were bracketed. Since the results for a = 90° 


and 150° gave a good indication of the shapes of the concentration distributions 
for optimum conditions, it could be concluded from the concentration distributions 
for a = 120° and 135° that optimum conditions had been reached or very nearly ne 
approximated even though | they were not bracketed. Tee oe ee. 
Dual Jet Injections.—In order to determine the improvement which could a 
be obtained in mixing by using multiple injections, experiments were conducted 
= two diametrically opposed 90° jets flush with the pipe wall. Angles of 
_ injection other than 90° and injections with three or more jets were not studied _ 
experimentally for reasons given later. The two 90° jets issued into uniform 
turbulent pipe flow, and the same momentum ratio was used for both jets 
in an attempt to have : sy mmetry with respect to » the horizontal diameter. 
‘The experiments with a single jet injection indicated that using the optimum = 
‘momentum ratio produced concentration distributions symmetrical about the 
pipe center line at large Z’s, so that optimum injection for a single jet was 
‘obtained when the momentum ratio was sufficient to carry the center of the 
injected tracer to the center line of the pipe. Thus, it was assumed that M,_ 
for the dual injection should be adjusted to produce jets which penetrate to 
the quarter points of the diameter. Based on this reasoning, on Eq. 10, and 
on M* = 0.013 for a single 90° Jet, runs were conducted with M, = 0.003 
and 0.004 for dual jets. However, the C, values were only slightly smaller 
7 - than would be expected for two wall sources. The apparent reason was that, 
_ in spite of the efforts even under laboratory conditions to achieve symmetry 
about the horizontal diameter, such symmetry was not achieved, and larger 
than expected mixing distances resulted. For M* > 0.006, the concentration — 
Ee indicated that optimum conditions for uniform flow were obtained | 
by causing the jets to penetrate somewhat past the quarter points of the diameter. 
With the maximum concentration between the quarter points an and tl the "center 


“center line. The reason for this behavior can be seen by visualizing a two-peak 
_ concentration distribution across the pipe diameter. Diffusion moves mass from 
4 the two regions of high concentration to the three regions of low concentration. © 
4 If the jets penetrate past the quarter points, diffusion from the large conerauatons q 
will ‘‘fill in’’ the distribution in the center of the pipe faster than at the walls” 
since there is less to “‘fill in.” Then diffusion will take place only toward the 
wall and the maximum concentration will be in the central portion of the _— 
Apparently, since it is not possible to obtain precise symmetry of the injection, 
: the practical optimum is for penetration past the quarter points of the diameter. a 
_ Thus, essentially the same concentration distributions and C, values were obtained — ; 
for the tests with M, from 0.006 to 0. 010, so that M* is in this range. For 
this optimum dual injection, Z,, = 0.28 and 0.17 for C, = 0.02 and 0. 10, 
respectively. The corresponding values of Z,, for two passive “aul sources 180° 
apart (9) would be 0.32 and 0.21, so that the use of the two jets gives only 


a 13% or a 19% reduction in Z,, (depending ont the definition used) as compared 


= 


J JET INJECTIONS 


‘indicate that the n mixing distance is is more : sensitive toa a source’s being off fll 
center a short distance than to being moved inward from the pipe wall a short — 
distance if a small C, (e.g., 0.01 or 0.02) is used in defining Z,,. This analytical 
_ indication of the effect which a slight misplacement of a center line point source 
_ would have on the mixing distance for small C, is an indication that the mixing 
distance for a jet with M, near M* would also be rather sensitive to the 
concentration distribution immediately downstream from the initial mixing zone 
and therefore to the transverse alinement of the jet. Sensitivity tests per ad 
™ ere not run for the jet-type injection, but during the experimental program, 
‘it became apparent that the mixing was sensitive to jet alinement. Even with 
“the utmost attention given to the alinement, it was difficult to obtain concentration - 
re symmetrical ; about the vertical diameter. With care, 
~ jets could be alined well enough so that the slight asymmetry was not 
evident in the C, versus Z graphs for C, > 0.03 for the single jet injections. | 
The distance which a concentration distribution can deviate from the center 
fer the asymmetry effect becomes significant could not be determined 
ete without measuring the concentration distributions i in more detail than — 3 
was done for this : study. However, the “results for an optimum 90° ‘Single jet yr 
in the present work were compared to those of an optimum 90° single. jet which 
was Slightly off center (7) and which, for optimum conditions, produced a po 
concentration located at r/D, = 0.05 on the horizontal diameter at z/D, = 
4, The optimum Z,, for C, = 0.03 was 55% larger for the previous study than 7 
for the optimum injection in this study. Also, was 17% larger 
than the value of 0 0.013 obtained in this study. 
7 ‘Tests with Secondary Currents in Ambient Flew. —The symmetry v which has 
been considered requires symmetry of both the flow and the injection. Since 
_ secondary currents exist in many pipes, the flow may not be symmetrical with 
to the center line and even a high degree of in alining the 


mixing as measured in the Thus, the possible of 
. currents on the mixing was considered further by placing a three-blade, fixed, — 
jection motor propeller in the pipe 13.5 diameters upstream of the point a. 
injection. The resulting secondary current consisted of a single cell swirl Stich 
was approximately symmetrical with respect to the e center | line and v which caused — 
a full rotation of the flow in approx 8diameters.§ == 
_ A single 90° jet injection (M, = 0.013, Run G1) and dual 90° jet injections 
sg = 0.010, Run G2) were used to compare the mixing characteristics with | 
and without the secondary current. The semilog graph of C, versus z/D, is 
"presented in Fig. 5. Variable C, was plotted versus z/D, rather than Z because 
“the propeller changed the structure of the pipe flow and thus probably changed 
_e, compared to the previous tests. Experiments \ were not conducted to determine 
A comparison (Fig. 5) of the single 90° jet with and without the secondary 7 
current shows that the swirling motion caused the log C, versus z/D, curve | 
to break away from the ‘‘symmetrical’’ slope sooner than for the condition a 
with no swirl. The concentration distributions showed that the swirl apparently © if 


two passive sources for uniform flow. 
- _ Effects of Jet Misalinement.—Fig. 1 shows C, calculated from the analytical : 
q | 
| 
| 


distribution. Thus, even though the propeller probably increased the 
increased the mixing distance. On the other hand, a comparison of the C, values 
for two 90° jets (M, = 0.010) issuing into a flow with and without the secondary 
current (Fig. 5) indicated that the swirl decreased the mixing distance for - 
dual jet arrangement. There is too little data to obtain a definite explanation = 
_ for this observed decrease in mixing distance, but the explanation might be : 
as follows: With no secondary currents, the two jets approached each other _ 
and the two concentration distributions which resulted from the initial mixing 4 
merged together about the center line. With the single cell swirl in the flow, Fi 


the top and bottom jets were deflected in opposite directions so that — did 
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Pest FIG. 5.—Costticient of Variation with Secondary Current in Pipe Flow a 
not merge at the center and therefore ‘produced a greater amount of initial 
| at i.e. a more ‘‘spread out” concentration distribution at the end of the 
_ Uniform Flow.—The results from the dual 90° jets did not indicate a significant _ 
_ decrease in the mixing distance as compared to the single 90° jet for uniform 
_ flow. It was assumed that any decrease which could be obtained by using two > 
_ jets with a > 90° instead of one jet would be on the same order of magnitude - 
_ as the decrease for a = 90. Therefore, dual jet injections for a > 90° did 


| 
= 


thew is only a small reduction in the mixing distance, 


3 using three sources instead of two, and that there is no reduction in Z,, when 
_ ~using four or more wall sources instead of three. In view of these results, 
andi in view of the fact that: (1) The optimum condition for two jets corresponded 
to the concentration distributions’ merging at the center of the pipe; and (2) 
_ the optimum two-jet injection gave only a small reduction in Z,, compared = 
to o the optimum single jet, it appeared that there would probably be 


currents jellies than being truly unifor orm. Because of the various types of an 
7 flows which can exist, it is difficult to generalize on possible interactions of 
jet injections and secondary currents. Nevertheless, some secondary currents 


for the optimum conditions with uniform flow. Even though | the us use of three _ 
or four injections apparently would give only small reductions in Z,,, compared 
to one or two injections for uniform flow, Fig. 5 shows that the use of two 
“injections provides a definite advantage compared to one injection in a pesticnion 

flow with a secondary current. At least for some flows, it might be expected © 

- that th the use of two or more injections would provide a sort of factor of safety 

against possible. effects of secondary currents on Zn there is not 
> 


Pump Mixing. mixing deserves some even this 
has not dealt with that subject, and even though there is apparently very little 
available data. Clayton (2) has shown that nearly complete mixing (less than 
a variation in concentration) was obtained in a mixed-flow pump with a 20-in. 7 
discharge line and with a flow of 8.2 cfs. The injection wa: was made at the bell- reaail O 
suction port, which was submerged in a sump. 
_ The efficiency of pumps gives one indication that p pumps should serve 
_ rapid-mixing chambers. If a pump has an efficiency of 80%, then 80% of + 


iz energy input goes into the mean flow while 20% is lost, with a significant 
part of the loss being in the generation of turbulence which produces the mixing. 
_ and, therefore, little mixing generated by the blades of a large turbine. _ 
a If a pump is being considered for mixing, the considerations should el 
- possible volitalization of the additive due to low pressures on the suction side — 
of pump and any possible corrosion of metal parts, deterioration of seals, etc., 
due to the concentration of the additive inthe pump. = = = 
7 yy Pipe Mixing.— jai —If speed of mixing is not ot important and if sufficient pipe length — 
_ is available, then a simple injection at or near the pipe wall may be all that — 
ts required. A distance corresponding to Z,, = 1.0 will produce enough mixing 
_ to give C, = 0.01 when a negligible density difference exists between the fluid 
Ly in the pipe and the injected fluid. The corresponding maximum deviation from 7 
the average concentration would be less than 2% (9). For a friction factor of 
002, A of 1.0 corresponds to 147 pipe diameters with Sc, = 1.0, or f = _ 
0. 01 gives 208 diameters for ; an 0. The use of two points sources, 180° 4 
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apart at the wall, would reduce Z,, from 1.0 to 0.37 for C, = 0.01 (9). a 
there is insufficient pipe length to accomplish the mixing by ambient turbulence, 
or if more rapid mixing is required, then jet injections could be considered. ; 
When it is desirable to use a combination of initial mixing due to the injection — ‘ 
and ambient mixing due to the pipe flow, the selection of the best injection 
ae ‘system for promoting rapid mixing in a specific application requires consideration 
of at least six factors: (1) The required degree of uniformity or the required — 
coefficient of variation for the additive; (2) the speed with which mixing needs _ 
_ to be accomplished; (3) the pipe length available for mixing; (4) the hydraulics — 
of the pipe flow, i.e., the friction factor and secondary current pattern; (5) 
the accessibility of the injection region and possible problems installing and = 


operating an injection system; and (6) the cost and availability of power for 


4 In general, the required degree and speed of mixing will depend on the particular 
7 agyeeaien. When using the method of dilution for measuring discharge in a 
- pipe, e.g., it may be advantageous to measure the tracer at only one point a 
in the pipe flow at a cross section where C, is small, say less than 0.01. For 
‘other applications, larger C, values 1 may be acceptable. 
_ The pipe length available for mixing should be determined. Some injection 
configurations may be excluded from consideration because the flow distance a: 
required to reach the specified C, is longer than the available ) pe length. For — 
established, uniform flows with negligible secondary currents, the results from — 
Figs. 3 and 5 may be used to estimate the distance required to obtain a acertain 
degree of mixing for the optimum momentum ratios for injections with one 
and two jets. For two jets, e.g., with M, = 0.010, C, = 0.10 would be achieved 
at Z = 0.17 and C, = 0.01 at Z = 0.33. For f = 0.015 with Sc, = 1.0, Z 
= 0.17 corresponds to z/D, = 29. Similarly, Z = 0.33 gives z/D, = 56. For z 
a 3-ft (0.91l-m) diam pipe with a — of 5 fps (1.52 m/s) and with 1-in. 
54-cm) diam injection tubes with M, = velocity and discharge 


= 1) 


= V0.010 


Power requirements are considered later. on 
Asn mentioned previously, many hove significant secondary currents 
which « can destroy the symmetry that is necessary to achieve the optimum rates 
of mixing obtained for uniform flows. Some experimental results were obtained — 
for a flow with a particular type of secondary current, but that current was 
a rather strong, single cell swirl and the results for that case are not necessarily 
i indicative of the behavior in other situations. On the basis of present | information, 


appear to be advisable to use : two ‘or more 90° jets ; equally “spaced around ~~ 
pipe periphery, w with each jet having M, of 0.006 to 0.010. The use of more 
than two jets provides a factor of safety against the unknown influences of 


| 


in . which P = power ‘in kinetic energy, of pipe flow; Q= discharge; V = = 


a velocity; D = diameter; p = pipe; j = jet; and N = number of jets. ‘The 
e variable P, eocmnin as D, increases for a given pipe flow and given . M,. Ea. 


additional p power requirement could come from the need to | pump the jet against a 
the internal pressure in high pressure lines. However, most of this additional 
power requirement can be eliminated by a closed system where the fluid for 
the jet injection is extracted from the pipe. Then a small-power, high-pressure 
_ pump could be used to inject a concentrated tracer (additive) into the closed : 
sy stem. ‘The ‘injection should be made in the by- pass line” ‘So that complete 
occurs before the jet enters the main pipe. This mixing could be 
DD samy by injecting upstream of the pump or by providing sufficient — 
tion of the by-pass line. Friction losses and po? efficiency must also be 
_ included i in a calculation of total power 
_Assume that two jets with 1- -in. (25. 4-mm) diam are to be used to obtain 
M, = - 0. 010 in a 10-ft (3.05-m) diam pipe with a flow velocity of 5 fps (1.52 
m/s). The value of P, would be 1690 ft-lb/sec or 3.1 hp (2.3 kw) for the 
jets. Using D, of 2 in. (50. 8 mm) would reduce the required power by a factor = 
_ of 2 for the same M, (Eq. 14). ‘eine: : 
- When injecting a miscible fluid into a flow in a pipe, the use of jet injections © 
can speed the mixing relative to that which would take place just due to the ; 
ambient flow. The jets c can be at the pipe wall and therefore do not require — 
any appurtenances « or ‘devices be placed inside the pipe. 
4 In this study, laboratory experiments were conducted to investigate the 
ei of jet-induced mixing in pipe flows. The coefficient of | variation, 
C,, of the concentration distribution was used as a measure of the degree 
4 of mixing at a given cross section. The dimensionless mixing distance, oe 
- was defined | as the flow ‘distance required to to ) obtain a a specified il Increases 
m i: the effectiveness of the jet mixing produced decreases in Z,, and | vice versa. an 
For uniform pipe flow with a given number and orientation of jets, the ratio, cf 
M,, of the jet momentum flux to the pipe-flow momentum flux is the parameter — 
indicative of the mixing induced by the jets. For a single jet oriented at various — 
- angles to the ambient flow, the momentum ratios for optimum mixing and the — 
corresponding m mixing distances were determined. Increasing the angle of injection “ 


- - _ secondary currents which might exist in the flow. For pipe flows with natural _ 
= | secondary currents, there is a sparsity of data on the behavior of jet injections 
and more research is needed on thistopic. 
Reg pment nowe we Pmen ate no the ne 
| 
q 


90° (cross flow) to 150° | (almost. counterflow) ‘decreased the oj optimum Bas 
but at the cost of a large increase in the momentum ratio and in the power 
requirement. Under many circumstances, it might not be practical to consider — 
any angle of injection other than cups: ony 
were also conducted with two diametrically opposed 90° jets 
in uniform flow. Optimum mixing was obtained with the momentum ratio for 
each jet being on the order of one-half of the optimum momentum ratio for 
= single jet injection. Thus, even though the use of two jets in uniform flow 
did not produce a significant decrease in the mixing distance compared to : 
single jet, the two jets would require less total power than a single jet. Furthermore 
limited experiments in a pipe flow with an induced secondary current indicated 
that the use of two or more jets for the injection would provide a factor of 
against the effects which could have on 
e work upon which this publication is based was supported in part by 
funds provided by the Office of Water Research and Technology, (Project 
A-091-111), United States Department of the Interior, Washington, D.C., 
_ authorized by the Water | Research a and ‘Development Act of 1978. Contents 
_ of this publication do ‘not necessarily reflect the views and policies »s of the Office 7 
of Water Research and Technology, United States Department of the Interior a 
nor does mention of trade names or commercial products constitute their 
"endorsement or recommendation for use by the United States oman il 
_ The research was conducted while the writers were affiliated with the Department 
of Civil Engineering of the University of Illinois at Urbana-Cham mpaign. 7 The 


project was through the University’s Resources Center. 
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INTEGRAL EQUATION 5 


Bruce Hunt’ and Lewis’ T. T. Isaacs’ 
Integral equation techniques provide one of the most accurate and efficient — = 
methods known for obtaining numerical solutions for certain types of harmonic © 
7 _ boundary-value problems. One of the biggest reasons for this relatively great . 
je - accuracy and efficiency is that the solution of the Laplace equation and its 
— awe conditions is reduced to the solution of a single integral equation — 
a i which the unknown is a function of one fewer variables than the unknown 
7 _in the original boundary-value problem. For example, a two-dimensional problem 
is reduced to the solution of a one-dimensional integral equation. Thus, the _ 
_ resulting set of simultaneous equations that must be solved contains N unknowns, 
_ whereas the solution of the problem with either finite « difference or finite element _ 
methods would require the solution of a set of enema, equations with 
approximately N? unknowns. 
Linear, Fredholm integral equations a are classified as the 
“and second kin kind, respectively 


The functions, f(x) and K(x, &), in Eqs. 1 and 2 are known, and the unknown 
function, occurs both outside and within the integral for a Fredholm 
‘integral equation of the second kind. Both types of integral equations are ss 
; encountered in applied problems, but equations of the second kind are usually 
preferred whenever a choice is available. This is partly because the theory 
bg for equations of the second kind is relatively complete, and mathematicians — 
have been able to ) apply the Fredholm theory to equations of the second kind Fi 
"Reader in Civ. . Engrg., Univ. of Canterbury, Christchurch, 1 New Zealand. 
oe? Sr. Lect. in Civ. Engrg., Univ. of Queensland, Brisbane, Australia, se 7 
Note. —Discussion open until March 1, 1982. To extend the closing date one month, Ae 
a written request must be filed with the Manager of Technical and Professional Publications, ie 
} ASCE. Manuscript was submitted for review for possible publication on November i} 7 
1980. This paper is part of the Journal of the Hydraulics Division, Proceedings of the 
American Society of Civil Engineers, © ASCE, Vol. 107, No. HY10, October, 1981. 
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a prove existence theorems for solutions to the Dirichlet and Neumann problems al 


on [Kellogg (2), and Sternberg and Smith (10)] . Probably ‘more important from 
ri engineer’s v iewpoint, though, is the fact that equations of the second kind — 


often provide more accurate and stable numerical solutions. This is because — 
_ the numerical solution of Eqs. 1 and 2 always reduces to the solution of a 
set of simultaneous, algebraic equations, and the presence of (x) outside the 
_ integral in Eq. 2 strengthens the main diagonal of the coefficient matrix. The 
net result is that equations of the second kind often yield a set of simultaneous 
equations that is well conditioned, and decreasing the nodal spacing (which 
_ increases the number of nodes and equations) leads to a more accurate solution - 
for (x). On the other hand, decreasing the nodal spacing in equations of the = 
first kind often leads to a poorly conditioned set of ‘equations with a resulting 
i in accuracy in the solution for (x). 
_ Techniques for the use of integral equations in irrotational flow problems 


have been investigated fairly intensively since the well-known work by Smith 
and Pierce (9) in 1958. On the other hand, integral equation methods have _ 
~~ been applied to ground- water problems since the work by Liggett (6) = 
1977. This is largely a result of the fact that irrotational flow problems can 
be formulated either @ as s Neumann or Dirichlet problems, both of which lead 


similar methods to ‘the mixed boundary- value problems in . ground- water flow 
often leads to a mixed system of integral equations of the first and second 
kinds. For example, Liggett (6) used the following equation in his work when 
‘. (x,y) wasa point on the boundary where the boundary tangent turned continuously . 
_ In Eq. 3, @ = velocity potential; s = arc length along the ‘boundary; r= distance 
between the “fixed” point, (x, and the integration point (é,n); and n 
arc length in the direction of the outward normal. When (x,y) is a point on 
7 a boundary of specified, a¢/an, Eq. 3 becomes an integral equation of the 
_ second kind. However, an integral equation of the first kind is obtained when 
(x,y) becomes a point on a boundary of specified o. The implication is that 
the resulting set of algebraic equations that are used to obtain an approximate, 
numerical solution of Eq. 3 may become ill- conditioned when nodal spacings — 
become too close on boundaries of specified . In fact, Liggett (6) mentions | 
that he has encountered an ill-conditioned matrix when nodes are spaced too — 
closely near points where and a¢/an change rapidly. Later work by Liggett 


ds 


_ The purpose of this work is to investigate an alternative integral equation 
formulation for the two-dimensional, mixed problems of ground-water flow. _ 
This formulation yields a simultaneous set of integral equations of the second 
kind which hhave boundary values of the velocity potential and the stream n function 
& its unknowns. The formulation has advantages of simplicity and compactness, ; 
and it is hoped that the set of integral equations will yield a set of algebraic 
- equations for its approximate solution that remains well- conditioned when 

boundary node are decreased. 
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INTEGRAL EQUATION FORMULATION 


Let o(x, y) and v(x, be and stream a 
tyne ee flow problem. The exterior boundary of the flow region is denoted A 


"problem i is equivalent to finding a potential 

ia which V=1; 2 genet 


iy; and w(z) = a function 


= analytic within I. It is idamadl that either or \ are specified on the adjoining ~ 


FIG. for Solution of Bound 
Cauchy integral theorem ICarrier, 1 Pearson (1)] that 
4 


= coordinate of a point within ¢ = coordinate of a point on 
a and the contour integral is calculated around the closed boundary I by traversing ¥ 
4 _ so that the flow region lies on the left. A definition sketch for this problem — 
Se... integral equation i: is obtained from Eq. 5 Y= letting z approach a point 
te the boundary. Carrier, Krook, and Pearson (1) show that this can be done ~ 
* placing z on the boundary and then deforming the boundary near z into < 
a portion of a small circular arc that places z within the flow region, as shown es 


in 2. end letting the radius of the circular arc 


pa 


in which (P) denotes the Cauchy principal value of the integral; -and 2 = = interior 
angle between the boundary tangents at point, z. At points where the boundary 7 
Eq. 6 can be rewritten as two real e equations by Setting w = + ib and 
by taking the logarithmic derivative of z=r in which r = 


la 


vwollo\sds om’ 


5 

- Eq. 7 becomes an integral equation of the second kind when (x, y) is a point ; 
on a boundary of specified, ¥, and Eq. 8 becomes an equation of the second 
kind when (x,y) is a point on a boundary of specified o. Thus, the solution 
of Eq. 6 for a mixed problem is equivalent to the simultaneous solution of 
_ two sets of real, integral equations of the second kind with 4 and as unknowns. 
Eq. 6, rather than Eqs. 7-8, will be used herein since the use of complex 


variables makes a considerable ‘in both the mathematics and 
coming. ate 4 


| 
| 
j 
| 
and @ = argument of 1 — 


n INTEGRAL EQUATION FORMULATION 1201 

The numerical solution | of Eq. 6 proceeds by first dividing the boundary, 

I, into N segments. These segments should be shorter in regions where 

or & change rapidly and longer in regions of slowly varying or . It is convenient 

_ to assume that w(t) varies as a linear function of the complex variable, t 
along each segment in a piecewise continuous manner. Thus, along the segment 
PUP US 


joining the nodes, , w(t) is approximated v with 
9 can be inserted into Eq. 


+ 
dt = 1+ in 
\ 


by vsing ¥ to 


ultiple- viees, and care must be taken to ensure 


‘The logarithms i in in Eq. 10: are 


3.—Beundery Geometry for integration Along Adjacent t to Singular 
Node wh <= 7; and px 


fot eh Te 
that a correct branch is used for all calculations. eye the imaginary part of re) 


logarithms is just the change in of t varies from t- 


Eq. 10 cannot be used to compute ‘the that comes from 


- over the elements that lie on each side of the singular node at z = - 
this. case, calculation o of the Cauchy value gives 


dt = Wiar “1 + + i(a — 8) 


3, = — between the boundary 


> 
7 
7 


on lines joining points, ¢,, 


baits and q 

The angle, @, can be computed numerically in a computer by first calculating fl 3 
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0= 


in p< if B > 0; 6 = — B; if B < 0, 

The difference, — 6, vanishes if the boundary elements on each side 

of the node, z = ¢,, are ‘straight- line Segments. 


- It is important to point out that Eqs. 10, 11, 


12 hold for boundary elements — 
that have any geometry. The fact that the integrals depend only upon quantities 


calculated at the endpoints ; of each element is a direct consequence of a theorem : 
in complex variable theory that states that the value of a contour integral of 
an analytic function depends upon the beginning and | ending point and not upon 
path followed between these two points. 
Anisotropic 


The seepage head, A, is defined by by 


any of water; and Fn cceleration due to 


If the principal axes of seepage are parallel to the x- and y-axes and ‘. 


_k, are the permeabilities for directions parallel to the x- and y-axes respectively, 


anisotropic problem may be converted problem on the 


~The p perme for the transformed plane i is defined by 


They oleae potential function, , is defined by 


and iil stream function is defined, in the normal manner 
. The preceding computational algorithm for isotropic seepage is “used in the 


transformed plane to solve the anisotropic problem. 

An advantage of the transformation used is that the solution method can 

be easily extended to multizone, anisotropic 7 


| 


a 
| 
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Mutrzone, Anisotropic Seerace Yo, 9 gringo .bodinn 
_ Sometimes the region of analysis consists of a number of zones such that 
k,, k, are constant within a zone but the values vary from one zone to the — 
next. "The _ anisotropic seepage within any ‘one zone can be analyzed by the 
_ method given in the preceding section but the overall solution needs compatibility 
_ At a point on a boundary between zones a and b, fA has the csi 


—. 


a b (20) 


The numerical solution of Eq. 6 carried out using N to 
the boundary into discrete segments. Then use of Eq. 11 to calculate the ; ; 
yong from the two elements adjacent to the singular node at z=t, 
ps owt lige 9 od odT 


in which the integrals are computed f from Eq. 10; and 
= t,. The angle, 2, no longer appears in 21 because of cancellations 
occur when Eq. is substituted intoEq.6. od} 
_ Eq. 21 may be solved in two different ways. The first method uses a Gauss-Seidel a 
"4 - iterative procedure in which the latest available estimates for w are substituted 
4 into the right side of Eq. . 21 to obtain the next estimate for w,. Then the previous 
“estimate for , is replaced with the real part of w, on boundaries of specified, 7 
: , and the previous estimate for \, is replaced with the imaginary part of w, _ 
_ on boundaries of specified, . The values of | on each streamline boundary ia 
are set equal to the imaginary part of w calculated at the first node on that 
Streamline for each cycle. The modulus of the difference between values of 
w calculated for two successive cycles at each node is summed for each cycle, — 
and calculations are stopped when this sum ceases es to change significantly from 


The iterative procedure can be ote. coded into a computer Program which 
ea 0 satisfy Eqs. 19 ‘and 20 before proceeding to the next zon zone. e. The = 
fact that the iterative method converges in most cases indicates the good 


conditioning of the equations. However, problems occur with very high ratios 
of k, for which the ‘second method of ‘solution i is found to » be satisfactory. < 
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‘The ‘second method requires assembly of all equations for all zones 


_inwhich G= 8, + ab 
A=a, + ia, = G|1+{— 
and =6,+ib,=G ——— }l | —— 
_ Eqs. 24 and 25 are derived te = pai for the integration along a segment 
‘y The complex equation, Eq. 22, may be split into two equations, one each 
for the real and imaginary parts. Multizone problems are formulated with h 
as the nodal variables instead of and dy because may vary from 
_ 


zone to zone. The real equation is 


Dd + + bk + = 0 
and the i imaginary equation od wt pl 19200, 


_ When a node lies on a beuadeny common to two or more zones, the fi Final 
equations for A and w at that node are obtained by combination of the oe 


equations from each zone. for al wi 
‘The full set of e 


inv which [C] = a matrix composed of known coefficients; and {a} = a vector _ 
comprising the nodal value of and arranged so that a,, m and 
Matrix [C] is formed by accumulating the contributions from each ch segment 
and for all zones. The procedure can be implemented directly from Eqs. 26 
and 27. The coefficients for the real equation for node, i, are stored in row 


- lof and the coefficients for the imaginary equation are stored in 


| 
q 
4 


two segments adjacent to node, i, and the summations represent the contributions 
_ from all other segments around the zone boundary. The column position for 


. each coefficient is determined from the location of the variable in {a}, e.g., 


considering the term, a,$,,,, in Eq. 26, a, is added into location [2i — eo 
- 2¢j7 + 1] in [C]. This method for the assembly of equations is analogous — 
i to the assembly of equation in the finite element method. In the analogy, each 
zone is treated as a separate element. ofT 
a After assembly, the full set of equations is modified to take account of the © 
| - boundary conditions. If h,(or ),) is known, the equation for h, (or ,) is replaced 7 


by an equation which specifies 1. A, (or 1. §,) = known value. If i is a point 
on a boundary where | = q but q is unknown, one nodal value of J on that 

boundary, e.g., v,. is chosen as the unknown and the equations for the remaining _ ; 
values of are replaced by 1., — 1.,, = 


aru 


EXAMPLES 


_ from the analytical solution for a semi-infinite region so that the results may . 

be compared with analytical results. The cutoff wall is of zero thickness and, — 

as a result, if only one zone were used for the analysis, numerical problems _ 
would arise in the calculation of coefficients when z is located at a node on | 


the wall. These problems are avoided if the region is subdivided into two zones 
_ as shown in Fig. 4. Since w(t) is assumed to be a linear function of t, a reduced 


_ The variation of h along the dam and cutoff wall is shown in Fig. 5, where 
Us _ the calculated and analytical values are compared. Lafe, Montes, Cheng, Liggett, a 
= and Liu (3) use a similar example to show how specialized quadrature formulae _ 


i, =: row 2i. The first five terms in Eqs. 26 and 27 are the contributions from the ar : 
_ -The first example (see Fig. 4) is of flow beneath a dam with a sheet-pile, 
cutoff wall. The lower boundary has been chosen to coincide with a streamline 
m3 
2 
a _ spacing of nodes is employed near the singularities to improve the —_§# 
P 
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encountered at the cutoff wall tip in this example. Undoubtedly, a 
7 would improve the accuracy of this solution. However, the agreement between 
=e and calculated solutions in Fig. 5 is close enough to suggest that an a 


accuracy sufficient for most engineering purposes can be obtained with no special 


_ treatment other than closer node spacings near singularities. The calculated 


: value of g/ KH is 0.464 and has an error of 1% when compared with the analytical 


—Seepage Heads along Sheet-Pile Wall and Dam Base in Example 1 — 


® 


ag 
ba 


An alternative solution is avai 


7 _ are shown in Fig. 7. Since the purpose of the original analysis was the determination 

of pore pressures for a stability analysis, a comparison of the calculated pore - 

7 pressures along AB and CD (see Figs. 6, 7) is presented in Fig. 8. The seepage _ 
discharge is also available from both analyses, and the values differ by 2%. 


This second example demonstrates the practicality of the method, and the | 


7 
| 
| 
‘ow The second example is the analysis of seepage through the mine tailings dam 7 
: _ shown in Fig. 6. The ratio of k,/k, is 4.0 in all three zones and the values 
fork, were k,,:k,,:k,, 1:4:30. 
> 
1EM, 
FIG. 6.—Example 2: Cross Section of Tailings Dam 
lable as the dam had previously been analyzed . 
. using the finite element method. A 6 node triangular element was used in the 
finite element solution and the mesh comprised 172 elements and 399 nodes. 
A total of 37 nodes is used in the integral equation solution. Of these, 5 are : 
7 located on the common boundary between zones | and 2 and 7 on the boundary 
between zones 2 and 3. The positions of the phreatic line from the two analyses > 


FIG. 7.—Location of Phreatic Line in Example 2 - BIEM, 


FIG. 8 —Pore Pressure Heads AB and CD in Example RA 


differences between the finite element solution and integral 


_ A major objective of the study was ‘the development of a computer program . 
which could be easily used in practical applications. Since the method can model 
_ irregular boundary geometries and a number of different anisotropic zones, 
_ this objective has been achieved. The second example demonstrates one s such 7 
application. Furthermore, the method has a ‘significant practical advantage o over 
. . ss finite element method because of the big reduction in data preparation. — 
- This advantage is still important when values of the solution are required within © 
_ the region of analysis because only the points where the solution is — 
need be specified. The value of w at any specified internal point can be evaluated © 
from Eq. 5 with the integral | taken around the boundary of the zone containing 
Unsteady, free-surface solutions have not yet been attempted with this” 
_ formulation, but it seems reasonable to expect that they could be carried out - 
ina a way. This is because an an | unsteady-fl flow Problem is solved 
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numerically by cz calculating sc solutions number c of boundary- 
BS problems in which the free surface position and shape is altered with 
time. ‘The position of the free ‘Surface for each problem is controlled by the 


lerivative 


values of along the free surface, it should be fairly easy to extend the technique 


L to obtain unsteady, free-surface solutions in two dimensions. On the other hand, 7 


three-dimensional and axisymmetric problems cannot be attacked with | the 
-variable approach considered herein. 


easily coded into a general purpose computer program for the analysis of — 
_ multizone, anisotropic, two-dimensional flows in porous media. This formation — 
leads to a well- conditioned Set of equations and allows nodes to be ddesdty 
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The following symbols a are used in this paper; 


™ of singularities. The program can be readily applied to practical problems. __ 7 | . 
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k ‘permeability; 
= length in direction of outward normal; 
pressure; 
_ distance between z and 1; ANCE OF 


arc length along boundary; sed € 
point (é, n) on boundary, I; 


components; 
+ iy, complex potential; 
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exterior boundary of flow ‘oscil Tt 


density of water; anes of 


velocity potential function; 
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_ Erosion of cohesive soils has been studied extensively from is hydraulics 
"point of view but less attention has been | given | so far to s soil behavior. The 
‘importance of understanding the soil strength mechanisms controlling erosion | 
“resistance i is recognized (27). Recent work has shown that rate process theory _ 
-may be useful for understanding the erosion resistance of cohesive soils and > 
emphasized the importance of physicochemical aspects such as pore and eroding 
The purpose of this | paper is to of designed to 
the applicability of both rate process and double layer theories to the erosion 
of cohesive soils. The analysis will show that critical shear stress for surface 
4 erosion of a cohesive soil exhibits the same dependency on the number of 
- interparticle bonds per unit area as reported for soil strength at higher stress a 


levels indicatin that they are phenomenologicall similar, 


_ There is an extensive body of literature dealing with the erodibility . of cohesive - 
soils and only a brief review of aspects touching on this study will be included. , 

7 To that end, only work dealing with weak, saturated cohesive soils will ill be 
The fact that a sediment exposed to fluid stressing is eroded i is ; well known 
P o The average fluid shear stress imposed on the soil surface can be related to is 
) the velocity of flow, and velocity and stress are often used interchangeably vi 
i in erosion studies. Since, according to Partheniades and Paaswell (25), the major _ 


> 


hydraulic parameter is the average fluid ‘Shear : stress, it is logical that relations — 
_ between the fluid shear stress that cause erosion and sediment shear strength 
would be sought. For weak soils, this approach has resulted in, at best, poor 


correlations, and Partheniades and Paaswell stated that there is no correlation. = 
Moreover, they suggested that the so-called critical shear stress that initiates @ 
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For cohesive soils, they concluded that surface erosion is controlled by 
oy interparticle bonds do actually control rates of erosion as Partheniades 
and Paaswell stated, then the concepts of soil deformation and strength based a 

x on interparticle bonds (14) may apply. Mitchell and his colleagues (12,13,14,15, 17) 
Me have developed a comprehensive theory of soil strength starting from r. rate process 
theory. For rate process theory to apply the rate to be described must B 4 i 
"temperature dependent. Several studies have shown that temperature does affect 
= rate of erosion of cohesive sediments (5,9,11,26). In fact, rate process teeny 
has been used by Raudkivi and Hutchison (26), Christensen and Das (5) iad 
’ Kelly, et al. (1D), to explain the temperature dependency of measured erosion < 
os rates. Christensen and Das and Kelly, et al., recognized the similarity of erosion 
_ at constant stress to creep at constant stress and used the methodology developed 
a: Mitchell, et al., to determine experimental activation energies and - 


caine and Das studied erosion rates by maintaining a constant flow 


of water through a brass tube tube lined with clay. _ They tested | kaolinitic and ite 
- clays and mixtures of these clays and Ottawa sand. They reported that both 
the stress and temperature dependency of erosion rates were in agreement with — -% 
rate process theory, although the experimental activation energies and flow a 
_volumes they determined were somewhat lower than those generally reported — 
for soil creep (14). Paaswell (22) commented on the potential usefulness of 
ae erosional model based on rate process theory, but cautioned that pipe flow 


_ Raudkivi and Hutchison (26) reported results of a flume study of kaolinite 
- erosion in which the effects of salinity and temperature were both evaluated. Pas 
- Although there is considerable scatter in their data, they concluded that, at 


no effect of temperature on erosion rates. erosion rates 
_ Reinterpretation of Raudkivi and Hutchison’s results indicates experimental 
- activation energies in the range of 5-10 Kcal / mole. These values are lower 
| than those reported by Raudkivi and Hutchison, but as they experienced _ 
considerable scatter in their data alternate interpretations are possible. Their 
_ plots of erosion rate versus temperature are U-shaped, which implies negative s 
f activation energies. It should be noted that Raudkivi and Hutchison did not of 
try to explain the mechanism of erosion resistance, but rather the observed 
Some attention must be given to the effect of temperature on the eroding 
wae Taylor and Vanoni (29) reported that an increase in temperature may — 
z. increase the intensity of the high-intensity turbulence fluctuations at the bed 4 
which accomplish the dislodgement and transport of bed particles. Vanoni and 
Brooks (30) reported that sediment discharge decreases with increasing water > 
temperature and that the concentration of suspended sediment decreases with © As 
increasing temperature. -Abou-seida and Arafa (2) reported that an increase in 
_ water temperature reduced bed resistance as well as bed form resistance. Both © 
the Vanoni and Brooks and that by Abou-seida and Arafa were for 


A in bed stress with ‘increased t temperature follows fro from = 


ef fect of temperature on viscosity and its influence on the shear stress coefficient . 
' An increase in water _ temperature of 10°K in the normal al range of water 
temperatures r results in approx a 25% reduction in viscosity. Since the ‘shear 
coefficient is approximately proportional to the inverse of viscosity to the negative | 
fifth power for turbulent flow, the result is about a 5% decrease in shear stress. { 
tT herefore, although changes in viscosity with temperature are significant, the a 
change in shear stress is much less and would tend to decrease erosion rates 
Physicochemical factors have also been down to be ae portant in erosion | 
_ (3,4,23,24,25,28). For soil behavior, Mitchell indicates that the primary role 
of the physicochemical forces is probably to control the initial soil fabric following _ 
sedimentation and to modify the normal interparticle forces” from what they 
would be due to the applied stresses. _ For surface erosion, it is reasonable — 


to assume that physico- -~chemical f orces. are the primary interparticle f forces since 


the a applied s stresses may y be, essentially, 


The rate theory (7) been applied to a variety 


particulate systems undergoing time-dependent flow or deformation. The theory 


P 


does not provide an exact description of the deformation mechanism, since - sa 
_ the basis of the theory is the assumption that atoms and molecules make up a 
flow units which are continuously attempting to move but are constrained from = 
_ movement by energy barriers separating equilibrium positions. The energy a 
flow unit must acquire to cross an energy barrier is assumed to come from 
thermal energy and applied directed ‘potentials. A ‘complete development 
the theory can be found in Glasstone, et al. (8). ter OF) 
a Applications of rate process theory to soil behavior are extensive and include 
_ work by Murayama and Shibata (18,19,20), Christensen and Wu (6), Mitchell 
(12, 14), Andersland and Douglas (1), Wu, Resindez, and Neukirchner Gl), ; 
Mitchell, Campanella, and Singh (15), Christensen and Das (5), Raudkivi .% 
(26), and Noble and Demirel (21). woils OOM 
a apply rate process theory to the erosion of a cohesive soil a deformation 
mechanism must be postulated. For the application of rate process theory to _ 
~ soil creep, Mitchell argues that the rate that flow units cross energy barriers — 
is to axial ‘Strain rate. For erosion ofa a Soil with ¢ constant fabric 


ane 
and tow (V,) can 1 be written 
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_ = particles leave the soil surface and hence the rate of erosion. This is essentially _ . . 
i the approach taken by Christensen and Das and is the approach taken in tis 
experimental =f 
g | 


ty 
(2 
_ These expressions al are re completely analogous: to those used in studies of soil 
creep, except that here the stress, +, is the average bed shear stress and the 
"rate is the erosion rate. In these expressions, , Ris the universal gas c constant, 
vin 
ki is Boltzman’s constant, and T is absolute temperature. From these | expressions 
‘it can be seen that it is the ratios of rates of erosion at different temperatures 
and stresses and not the absolute magnitudes that will be important in determining a 
_ Two approaches are possible for determining the tin'eh process parameters. 7 


First, tests could be run on identical cone at different ; temperatures and 


samples. second ‘method is to “increment or stresses and to > 
Sa experimental activation energies or flow volumes using equations (1) 
: and (2). In this study, rather than compute rates immediately before and after — 
an increment, the equilibrium rate at several stresses or temperatures is fitte 


-inaleast Squares sense "sing equations (1) and (2) to obtain experimental activation *~ 

heen. —Samples were tested ir in a water tunnel and designed and f; abricated 


33 ft (10 m) long. Water i is recirculated by a a 9-in. (23-cm) diameter four- -bladed, — 


dynamically- -controllable, ducted "propeller. The propeller is driven (through a 
- variable speed drive; 500-1,500 rpm) by a 10-hp electric motor. The tunnel cas 
_ is described in more detail by Gularte (10) 


Average velocities from 0-90 cm/s can be produced in the tunnel’s working © 


section, shown enlarged in Fig. 2. The velocity in the working section can 


be controlled by either the pitch or speed of the propeller or by adjusting 
a valve. In the working section, which has a “square 6-in X 6-in. (15.24-cm - 
x 15.24-cm) cross section about 10 ft (3 m) long, the entire sample is visible 
z= both sides through 2-in. (5-cm) thick transparent windows. The entire 
water tunnel is refrigerated and temperature in the tunnel can be maintained 
within +0.2° K for extended of or increased or decreased 
Two types of test were run; in one, temperature was held constant and velocity 
7 varied while, in the other, velocity was held constant and temperature varied. 
The tests were run on remolded Grundite (illite clay) at four different salinities 
(ranging from approx 2.5-10 ppt NaCl) and five different water contents (ranging © 
from 40-80%). In all cases, the eroding and interstitial fluids were maintained — 
at the same salinity and at a pH of 8.5. The remolded samples, which weighed _ 
a ‘egpeccinataly 44 kg with 1.5 sq ft (0.14 m’) exposed to fluid stressing, were — 


placed in clear lucite trays and inserted into the working section of the water = 


Samples were to an instrumented ‘shear The instrumented 


; ft (5.5 m) long, 6.5 ft (2 m) wide and 4.9 ft (1.5 m) high. Water, 38 gal (45 

- L) is recirculated in an essentially oval path (in the horizontal plane) about _ 
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asuring 
the average horizontal shear force on the eroding surface, from which the averngs 
; The center line velocity in the working section was measured continuously 
q with a ducted propeller-fiber optic current meter (10). The amount of material — 
im suspension at any time was measured with a calibrated laser- photocell 
arrangement which may be seenin Fig-22 
a _ Sample Preparation.—To prepare a sample, sodium chloride was added to 
metaphosphate) + was added and pH ages to 8. 5 by adding so sodium hydroxide 


: tap water to obtain water with a preselected salinity. Next, a dispersant (sodium 


FIG. 1 —Retrigerated Recirculating Water bra 


ana ‘AG. 2 Wetting Section of Water Tunnel 

pH was necessary to neutralize or - eliminate edge effects which would complicate 


“then added to the Grundite and mixed with a mechanical mixer to obtain samples 
Prior to placing samples in the support tray, they were remolded in their 

- storage containers. After placing in the tray and just before insertion into the 
tunnel, they were again remolded with a mechanical mixer. At this point, the 
= was considered to be completely remolded. Water contents (three) — 
were taken and undrained shear strength measured by a fall cone device just 
= placing and again after removing samples from the tunnel. After testing = 


| 
| 
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in the water tunnel, were stored in air-tight containers for 
Results from a typical tomperatere increment test in water are 
én in Fig. 3. Experimental activation energies were determined by least 
‘ squares regression | analyses of the data from each run. A summary plot is 
_ Shown i in Fig. 4. Results from | a ‘typical velocity increment water tunnel test 
are shown in Fig. 5. Experimental flow volumes were determined from the — 
slopes of semi-log plots of the rate of erosion versus bed shear stress, as shown 
in Fig. 6. Critical velocities were obtained by plotting erosion rates versus velocity 
7 and extrapolating back to the zero erosion rate (24). The. critical erosion stress 
was then calculated from the - critical velocity and the calibration curve for 
“the tray. Table | is a ‘tabulation of computed experimental activation | energies, 
flow volumes, and critical shear stresses with the corresponding test salinities 
ond 
_ The data in Table 1 may be compared with results from earlier studies on 
poll. Table 2 compares experimental activation energies determined from creep ; 
and viscosity tests, and erosion tests; the experimental activation energies for 
- erosion are somewhat lower than those for creep, but significantly higher than 
those for dilute clay suspensions. Table 3 compares experimental flow volumes — 
_ from this study and earlier studies. Flow volumes for Raudkivi and Hutchison 
and Partheniades were computed from their data and those for Christensen 
and Das by dividing their values” for bonds ‘per unit area by their assumed 
bond spacing. The basis for this “calculation may be found in Mitchell (14). 
Some significant trends are apparent in Tables 1 and 3. First, flow volumes 


soil creep and second, the tabulated experimental flow volumes for erosion 
An examination of the data in Table | indicates the | effects of water content 
and salinity on the rate process parameters. The variation of experimental flow 
—— with salinity exhibits a well-defined trend with experimental flow volumes 
decreasing roughly one-half over the range of salinities tested. The effect of 
_ water content alone is less clear; experimental flow volumes appear to decrease 
_ with increasing water content below the liquid limit (approx 55%) and then 
increase with increasing water content above the liquid limit. Experimental | 


activation energies exhibit no discernable trend with water content, but there — 


If the number of interparticle bonds per unit area is inversely proportional 
to the experimental flow volume as proposed by Mitchell then these results: 
“may be interpreted to mean that the number of interparticle bonds per unit 
_ area increases with increasing salinity and, above the liquid limit, with decreasing ~ : 
_ water content. Mitchell also showed how changes in the number of interparticle 7 
bonds per unit area for a normally consolidated clay could be accounted for 
by changes i in effective stress = increased effective stress meaning eed 
e the decreasing experimental flow volumes : 


with increasing salinity and eae: water content can be interpreted as 


= 
| appears to be a slight tendency for activation energies to increase with | 


4 q EROSION RESISTANCE 


FIG. 3. —Variation of Suspended "Concentration with Time and Temperature 
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FIG. 5.—Variation of Suspended ‘Sediment Concentration with Time and Velocity 
(1 in, = 2.t 54 cm) 
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clay particles with increasing salinity i is consistent with double-layer theory (id). 
The weak | dependence of flow volumes on water content is a reflection of 
the dominance of physicochemical effects at the soil water interface and the 
fact that only a relatively soft soil was tested. as i .on 

. Fig. 7 shows the variation of critical shear stress with water content. Lines 
are approximately parallel for the different salinities. The importance of salinity 
relative to water content is apparent. Fig. J is a a plot of critical shear stress 


4 ‘bonds per unit area a have been computed by dividing the assumed bond spacing © a 
@. 8 x 10~* cm) by the tabulated flow volumes (14). Also shown is the extension 


TABLE 3. —Experimental Flow Volumes _ 


= 


volume, | 
incubic 


wit 


(remolded) | erosion—water TS this study 


‘Tilite ar 
Kaolinite 

San Francisco Bay 
Kaolinite 


San | | 


mud 


(remolded) | j 


Illite (remolded) 
Sault St. Marie cay 


In Ref. 5. 
"In Ref. 
“In Ref. 6. 
Ref. 
“In Ref. 1. 


erosion—pipe 
erosion—pipe 


erosion—flume 


erosion— —flume 
ty 
creep 
creep 
creep 


ri 54-6.09 x 10°" 
: 
Christensen and Das* 
Christensen and ~ 
x Partheniades” 


ve 8. 10-"*| Mitchell, et al. 


A 


oft the line given b by Mitchell and Matsui (16) relating shear strength and number 
of interparticle bonds per unit area. The limits of the data presented by Mitchell | 


and Matsui are also shown approximately. It may be significant to note that — 


data given by Mitchell and Matsui for consolidated soils tend to fall along © 


line or above it given for clay suspensions tend to 


or below it. Aah 


_ The erosion results ‘anal here follow the same trends reported for ~ 
- higher stresses and clay pastes at high water contents. If the number 4 
4 interparticle bo bonds per ur unit area is controlled by the effective stress, ‘then the 


| Ses stress can be seen to be controlled in surface erosion by the interparticle 
Ses forces, which are described approximately by double- layer theory. 
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4 FIG. _8.—Critical Shear Stress and Shear Strength as a Function of the Number of 

> 
Alternate interpretations of these results : are possible. The : size of the e experi- 
s mental flow volumes could be used to argue that actual particles are the flow 
Units. However, the interpretation proposed here is reasonable and consistent _ 
with the data. More importantly, these results suggest that the concepts of | 
7 - soil strength a and def ormation proposed | to to explain s soil il strength and and creep behavior ay 


apply equally to surface erosion. pee 2770 iY . 2 


- Results from a series of tests designed to test the applicability of rate process - 
and theories for describing surface erosion of cohesive soils 


comparisons s with earlier s studies 


—_ 
& 
| 
&§ 
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: describing surface erosion of cohesive soils. The values of experimental activation 7 
energies for erosion are not significantly lower than those for soil creep. This ~ 
“ result leads to the conclusion that —_— bonds in both erosion and in 

} Flow volumes, a assumed to be i inoeieie related to the number of interparticle 5 
_ bonds per unit area, are shown to depend primarily on salinity. The resistance _ 
‘to erosion measured by the critical shear stress exhibits the same dependency me 
on the number of interparticle bonds per unit area as soil strength does. 
> take interpretation of surface erosion is consistent with concepts of soil strength. ' 
It also explains i in part why attempts to relate erosion resistance 
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CovER EFFECTS ON STREAM Flows 


By Y. ‘Lam Lau’ and Bommanna G. 


_ The presence of ice cover in a stream alters the flow characteristics to a _ 
oat extent. The normal flow depth increases due to the increased resistance _ 
resulting an additional ‘solid boundary and the average flow velocity 
decreases. In addition, both the velocity and shear stress distributions change. 7 
Therefore, one can expect differences in the momentum and mass transfer 

coefficients between ice-covered stream flows and open-water stream flows. _ = 


_ For a two-dimensional, fully developed turbulent flow in an open channel 


with a free surface, for which linear shear stress and logarithmic velocity — 
“distributions are often assumed to apply, the distributions of the kinematic 


_ viscosity, v, , and the mass transfer coefficient, Bus are parabolic, with a maximum 
at middepth and zero at the stream bed and at the free surface. MONTE EAS 
_ In analyzing the flow structure in an ice-covered stream, the flow is usually 


divided into two layers, separated by the line of maximum apse linear ; 


of maximum velocity. To avoid this obvious deficiency, the distribution of 4 
may be arbitrarily altered as by Shen and Harden (11) who adopted the eminetien 


that T, is constant in the central portion of the flow, as suggested by Ismail 


= 
= 


In addition to the arbitrary modification of the I, distribution, Shen 2 and 

7 Harden had to assume that the depth remains constant when the flow acquires 
an ice cover. This is also not correct because, given the same discharge, the 

“s appre will increase with the presence of an ice cover. In this paper, an alternate 2 


approach has been adopted to overcome the above difficulties. ~The ‘“‘k-e” 
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a 
ig 
q t, these distributions lead to parabolic distributions of v, and I, in —_ 
th a value of zero at the location of maximum velocity. This is obviously — 


to calculate the distribaticn, and v, distribution for two flows 
with the same given discharge, bed slope, and bottom roughness. One of these " rm 
y flows has a free surface at the top while the other has an ice cover of a given rs 
 .oughness. These results are then used in the two- dimensional mass transport 
equation to simulate the concentration distributions due to sources placed at 4 
different heights in the flow. The results give some indications of the effects 
of ice cover on the flow and on the vertical mixing. Sout a 


wy two » parameters, the energy of turbulent motion, 
k, and its rate of dissipation €. As a result, the turbulent er viscosity 
Eq. 1 permits the evaluation of the v, distribution from the distributions of — 
_k and €. The distributions of k and ¢, in turn, are determined by solving the 
- semiempirical transport equations for k and « along with the continuity and am 
_ For steady, two-dimensional channel flow, the equations of continuity, mo-— 
mentum and the transport equations for k and «€ take the following forms: 


au\ 


aud 


‘the channel bed and the y-axis is measured to the » x-axis in the 
_ vertical plane. Parameters u, v = the velocity components in the x and y directions, 
7 a “asap The parameter h = the flow depth; S = the slope of the channel 7 = 
Fe and C2 = constants; and G = turbulent 


q 


energy by the mean motion siven by: 


Liha 

_ in which $ = a scalar quantity such as temperature in the case of heat transfer 

and concentration in the case of mass transfer. The term 5, is is ‘the v ‘volumetric 

_ The empirical constants were determined by Launder oe Spalding (5) who — 
7 - _ considered a variety of flows such as flows in pipes, channels, mixing layers, 
: , jets ane wakes. The values arrived at by Launder and | Spalding are as follows: — 


(1) c, = 0.09; (2) o, = 1.00; (3) o, = 1.30; (4) c, = 1.43; 6) a= 1.92; 


and (6) 00 for mass transfer and 0.50 for heat transfer. 
The governing g equations listed in the foregoing are derived with the assumption 
_ that the flow is predominantly in one direction (namely, along x) and that the 
turbulent transports of u, k, €, and are negligible in that 
terms containing the second order derivatives in x are not included in the set 
of This latter assumption modifies the equations from elliptic type 


EITHER SURFACE OR ICE- te wt 


a 
to parobolic type. Parabolic equations are especially suitable for numerical : 
schemes which solve the equations for one station at a time and march forward 4 
along x. Such numerical schemes neater less computer memory and are executes 


¥ The fully developed seein of velocity and turbulent kinematic viscosity, - 
_— which are required for evaluating the effects of ice cover on the flow and : 
mixing characteristics, are obtained as asymptotic distributions resulting from T 
the solution of the governing equations. The same results could have been 
= by solving a set or ordinary differential equations representing the 7 
fully developed flows. Indeed, for fully developed | flows, the variation of flow - 
- : properties along x vanishes and therefore the governing equations become: Ne 


| 


dy 


| 
a 


_ The numerical schemes to solve the preceding system of ordinary differential 
equations are simpler in comparison to those schemes required to solve = 


partial differential equations. However, the prediction of concentration distribu- 
tions still | requires the solution of the partial differential equation (Eq. 6). In 7 

addition, the solution of the original set of partial differential equations can 

- be applied to the investigation of jet type discharges. Therefore, it was decided — 
solve the equations as expressed by Eqs. 2-7 the 


bal 


Es. 1-7 form a closed set of equations and, given values for | the empirical 
constants and for s,, therefore can be solved simultaneously for u, v, k, €, 
v,, G, and >. Boundary conditions for the above parameters have to be specified: 
(1) At the channel bed, ie., at y = 0 for all x; (2) at the upper boundary, 
i.€., at ys = h for all x =x and (3) at the initial cross section, i.e. .» at x=0 
for all y (initial profiles). Because of the parabolic nature of the governing 
equations, the boundary conditions at the downstream end of the calculation | 
domain need not be specified. The boundary conditions adopted | for the present 
calculations are similar to those used | by Rastogi and Rodi (8), | , except for the — 
—For the boundary layer type 
‘tas considered, “the: longitudinal velocity component near the channel bed, 
and near the ice cover universal law of the wall. 


and 
and the fl fluid densitypas 
« = the von-Karman constant, equal to 0.42; v = the kinematic viscosity of 
_ the fluid; E = a roughness parameter which takes a value of 9.0 for hydraulically — 
- smooth turbulent flows and (30.1/v, K,/v) for rough turbulent flows; and _- 
_= the distance from the solid boundary to the grid point nearest to the solid — 
boundary. The symbol K,, represents the size of the equivalent sand grain 
a The vertical velocity component v at solid boundaries is zero. | =e 
The turbulent kinetic energy, k, and its dissipation rate € near solid or a 
are evaluated using the following assumptions: (1) The } pyre of turbulent 
kinetic energy is equal to the rate of dissipation (i.e., G = €); (2) the turbulent 
shear stress t near the solid boundary is constant poh is equal to the boundary © i 


| 
| 
q 
= 


shear stress (i. e., T= t,); (3) the velocity distribution in the | vicinity of the 
solid boundary i is given by Eq. 9. 
_ Under these conditions, the kinetic energy of turbulence and its dissipation: 


rate solid boundaries can be evaluated as: 
‘yy 


For heat or mass transfer, or both, the flux a at solid boundaries is assumed 
> Boundary Conditions at Free Surface. —Following ‘the approach adopted by. 7 
Rastogi and Rodi (7), the free surface is treated as a symmetry plane for u, 
and Accordingly, the y gradients for u, k, and become zero. >». Rodi 
(8) has suggested that the presence of a free surface should have an hereclll 


of reducing the length ‘scale of turbulence near the surface and that a 1 symmetry 
boundary condition for € may not be completely satisfactory. Consequently, | 


it in , which k = the turbulent kinetic energy § all the free surface calculated using 7 
4 the symmetry boundary condition; and = the distance of the nearest 


point from the free surface in which k = k,. This condition is only | tentative 
. and needs further testing. However, the seulin show that it does not significantly 
affect the comparison between free surface and ice-covered flows. ied 
7 The vertical velocity component v at the free surface is zero. - {ide} 
oe For heat transfer problems the gradient of ¢ at the boundaries will be nonzero. 
However, this problem is not investigated here. 
ai 
Initial Profiles.—At the starting cross section, i.e., at x = 0, the distributions 
“of u, v, k, and € are not known a priori and, certain assumptions have to 
be made in specifying them. The distributions used in = present work are — 


} outlined in the following: 


rips 


2a in the vicinity of the channel bed, i.e. 9. In the remaining 
7 part of the flow region a uniform distribution is assumed. The assumed distribution 
has to satisfy the requirement that it yields the specified flow rate per unit 
width of the channel. The distributions for k and « are assumed to be uniform 
over the whole height of the flow field, having the values k,, and € » Tespectively, — 


_as given by Eq. 11. The vertical velocity component v is ‘assumed to mM zero 
_ 2, Ice-covered flow case: The velocity component u follows a et mic 
distribution both near the channel bed and near the ice cover. In the — 
_ region of the flow, u is assumed to be uniform. Here again, the assumed distribution : 
= to yield the specified flow rate per unit width of channel. Linear profiles — 
| k and « are used with the values k,, and ¢,, evaluated both at channel 
and at the ice cover using Eq. 11. ‘The vertical velocity vis” 


— 
‘ 
: 
i 
a Similar tO thal at the soli — 


ihe? 


(6) is adopted. For the sake of completeness, some of the salient features of aa 
_ the numerical scheme are presented in this paper. The forms of Egs. 3, i 
5, and 6 are such that one Single numerical scheme can be used to solve all 
hs the equations. Indeed, Eqs. 3, 4, and 5 can be expressed in the form of Eq. 
“~~ As an example, Eq. 3 is identical t to Eq. 6) when u is equated to > and a 
gdh/dx) to s,. Therefore, in presenting the | details of the numerical 
scheme here, only Eq.6isconsidered. 
y _ Inthe scheme proposed by Patankar and Spalding the fin inite difference equations 
. of Eq. 6 are arrived at by integrating the differential equation term by term 
over small control volumes. Certain assumptions are made Lge the variation — 
of along x and y - directions. . Referring to Fig. a), I is the upstream grid a 


id biloe ba) on os velin 


= 


| 4 


od tafe 


2.—(a) Control Volume; (b) Assumed Profile head in y Direction Grid 


located ata | distance Ax from the upstream od line where the values o 
are to be determined. The parameters J = t, J, and J + 1 are ‘grid lines in 
the y direction. J 1/2 and J — 1/2 a are the mid points between J and 
_ + land Jand J — 1, respectively. The control volume over which the equation 
q integrated is shown by hash lines. The variation of in the y direction 
is assumed | to be linear between grid points as shown i in Fig. 2(b). The variation © 
of b along | x is assumed to be stepwise. The step | change occurs right at the | 
upstream grid line. Between grid points ¢ is uniform and equal to the value — 
; corresponding to the downstream grid. Such an assumption of ¢ along x gives _ 
rise to an implicit form of the finite difference equations, 
5 When each term of Eq. 6 is integrated over the control volume shown in 
bs 2(a), the following relation is obtained: 


|. Transport velocities (u,v) and the eddy viscosity (v,) are defined along 


| | 
| | | 


e grid line (J) so that the equation for is linear. Further, the 
‘transport velocity, u, is assumed to be constant with y (i.e., upstream and 
g _ downstream faces of the control volume). The approximations introduced by _ “a 
this linearization and this representation of u, while not introducing an error 4 
for this uniform flow case, should be evaluated further before such a scheme — 

_ is applied to nonuniform or r unsteady, or both, flows. ong? clio civ 
(P + 0), Ay Ly ip 
proflie: of o so grid has imme, 


ia vero and Paine Of. buomes 


1/2 


Tet 


‘The integral of the source term S, over the ostiieil volume is ee. bd 


of Ss depend on the entity which ‘In Eq. 
_ allthe terms except the ones with s, and ¢ , represent the values at the downstream e 
- grid line J + 1. The subscript (J + dD is omitted for ease of writing. When - 
— Eg. 13 is written out for all values of J, a series of finite difference gel 


| 
(146) 
q 
| 
' _ several standard methods to solve such a system of algebraic equations. In q 


= the present calculations simple successive- ve-substitetion formulae proposed by a 


Patankar and Spalding w were used. avi 
AppucaTion oF NumericaL SCHEME tte 
Ps In order to evaluate the effects of an ice cover on the flow properties and 
the vertical mixing characteristics in channel flows, it is necessary to solve 
_ the governing equations for both the free surface case and the ice cover case. - 
_ For these two flows to be equivalent, it is assumed that both channels carry — m_ 
_ the same flow rate, have the same bottom slope and the same bottom roughness. | 
The flow depths, of course, will be different. The flow depth in the ice-covered 
i. flow will be larger than that of the free-surface flow because of the increased ot 
4 resistance due to the ice cover. Since the flow depths are not known a priori, a 
a trial and error method is al ™ procedure ete in the present 
- is illustrated in the following. hy 


FOR Free- Flow 


are assumed to be specified. A value for the flow 4 selected. 
(Subscript 0 is used to denote free-surface flow properties and i is used for 


4 (a) 


an 
| 


ICE COVER EFFECTS: 
, this flow depth is maintained constant in the longitudinal 
: Using this flow depth, a numerical grid is laid in the x-y plane as shown in be 
_ Fig. 3(a). Therefore, the problem becomes one of finding out what the slope — 
has to be in order to have uniform flow at the given flow depth. y barsveo-aal bos 
a. From the initial profile for u,, the shear velocity, v,, and thus: the ‘ty 
- values for k,, and «, at x = 0 are evaluated. The v, distribution is obtained 
_ from Eq. 1. Eqs. 2 and 3 are then solved using the numerical scheme to obtain 
the u and v profiles at the downstream station x = Ax. The value of the slope _ 
S$ in Eq. 3 is evaluated using v, and h,.  Rammmning uniform flow. ~~" ar 
3 the velocity ‘distribution at x = Ax, values. for and are 


used with Eqs. 2 and 3 to solve for the velocity distributions at the next downstream : - 
Station. This procedure is continued until a certain downstream distance is reached 
where the profiles of u and v, no longer c change from grid line to grid line, 
i.e., where the profiles become fully developed. At this point, the vertical velocity _ 
component is zero and the value of slope, S,, computed becomes invariant 
with respect to x. Therefore, for the specified flow rate per unit width g and 
& bed roughness K,, , the free-surface flow with bed slope S, will flow at. 
a uniform depth of ho. This completes the pe pen for open water flow. 
The values of flow rate per unit width of channel and the channel bed roughness 


for the flow under ice cover are the same as for the free-surface flow calculated 
The roughness of the ice-cover surface has to be specified. Let it be K,. 
‘The bed roughness for the ice-covered channel is denoted by K,. and its value 
_ is equal to Ca _ A value for the flow depth h, is selected, and ‘the computation 
’ is carried out using a a procedure s similar to that for the free-surface flow. — 
‘a computational grid for this case is shown in Fig. 3(b). The value of the slope _ 
‘Ss appearing in the momentum equation is evaluated using the shear velocities — 
at the the channel | bed and at the ice cover, v,,andv,,,ie. 
= computations are carried out until a certain downstream distance is reached ¥ 
where the profiles of u and v, no longer change from grid line to grid ome 
-i.e., where the profiles become fully developed. At this stage, the vertical velocity - 
component v, is zero and the shear velocities v,, and v,, become invariant 
with respect to x. Consequently, the slope S, computed using Eq. 15 is also 
_ invariant with respect to x. For the specified flow rate, g, bed roughness, ie 
and ice-cover roughness a4 the ice-covered flow becomes uniform with depth 
& for the slope S,. Now, if this slope S, is not the same as the slope S, of ry 
the uniform free surface flow previously computed, then a different value for 
_ h, is selected and the calculations are repeated until S, coincides with S,. The _ 


flow depth h, corresponding to this slope then yields an ice-covered flow equivalent © 
to the free-surface flow. ice: ler toed 


= 


Resutrs AND ANALYSIS j 


Using the the last section, fully developed profiles of 

- uand v, were obtained for three different flow conditions in both free-surface : 
and ice- ~covered channels. The hydraulic parameters for the three flow conditions 

_ are shown schematically in Fig. 4 and are also listed in Table 1. In the first 


J 
7 flow condition, denoted as run no. lI, the value for K, for both channels was _ 
q 


both channels. The ice cover was considered to act as an hydraulically : smooth 
_ In run no. 2, the bed roughness was increased to 5 mm and the ice-cover 
Ss was i the same as in run no. A In run no 3, me bed roughness = 


larger flow depths to transport the same flow rate with same bed slope. The 
increase in flow depth varied from 15%-31%. Calculations have also been made 
_ for a much deeper flow depth which is closer to a natural flow. In that case, 
Acs the depth increased from 2.7 m- -3. 0 m when an ice cover was present, an increase =z 
of 11.1%. The shear velocities at the bed of i ice-covered flows were smaller 
than those for the free-surface flows. The latter result indicates that if the _ 
bed of the stream is movable, then the sediment transporting capacity would 
_be diminished by the presence of ice cover. Furthermore, the bed forms resulting | 
3 from the movement of sediment would also be different from the free-surface 
flows, _ which in turn would alter the roughness characteristics. Therefore, the 
"presence ‘of an ice cover can n produce significant changes to the hydraulic 
characteristics of natural stream flows. However, it should be noted that if 
_ the ice-covered flow is divided into an upper and a lower layer, the Darcy 
~ friction factor associated with the lower layer is actually larger Ges & the 
‘Surface flow, even ‘though the bed shear stress has ‘decreased. This is 


Profiles of u “for the three runs are shown in Fig. 5. Solid lines represent _ 
the ice-covered flows and the dotted lines denote the free-surface flows. In 
all three runs, the ice-covered flow velocities are lower than those of the — re 
_ free-surface flows. The u profiles in ice-covered flows resemble the ones measured 
Oy Larsen (3,4) with close to zero gradients over substantial portions of the 
= region. Such profiles mean that the division of the flow into two layers, 
as has been done by previous investigators, can be very subjective. Velocity 
profiles in ice-covered flows are sensitive to the roughness characteristics of 
~ both bed and ice-cover. ‘When the bed roughness increases in comparison to 
- ice-cover roughness, the position of maximum velocity moves closer to — 
the smoother surface. profile becomes symmetrical when the roughness 
Velocity profiles for the free-surface flows are close to being logarithmic 
4 for about 75% of the flow depth for all three runs, but deviate from the logarithmic 
_ distribution for the remaining 25% near the free surface. This result is in agreement — 
with a a recent study by Song and Yang (12) who put forward the hypothesis 
7 ‘that the flow in a wide = channel consists of three different regions: a - 


§ = to. 5 mm. The resulting flow depths, shear velocities and slopes are listed in 7 
ah and are alco chown in Fio 4 exnected the ice-covered flows reauired 
} 
| > 
| 


aminar ar sublayer region, very close to the channel bed; an inner turbulent region , 

in the middle; and an outer turbulent region near the free surface. Based on | 
the experimental work of Vanoni (14) and re (9), mo and Fang argue 
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FIG. 4.—Schematic FIG. 5.—Velocity y Distributions 


that the logarithmic velocity distribution is applicable only to the inner turbulent 
~ region. For the outer turbulent region, they propose a parabolic distribution. _ 
__ A parabolic velocity distribution in conjunction with a linear distribution for 


= stress would give Tise to a constant | turbulent kinematic a when 


Z — 
ameters for Simulated Runs 
Run Number 2 Run Number3 
7 water covered 
22,240 22,220 
4 
- | 
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_ Gives rise to a parabolic distribution for v, as Shown earlier. Some experiments — 
carried out to measure the turbulent kinematic viscosity [Jobson (2)] do suggest _ 
“that the y, distribution deviates from the parabolic distribution, especially » il, 
‘the free surface which supports the suggestion that the velocity distribution © ( 
sod near the free surface indeed is not logarithmic. If it is not logarithmic, then o 
should it be parabolic as Song and Yang suggest? The answer to this question 
is not quite clear. Indeed, the v, distribution which is calculated from the k 
_ and ¢€ distributions using Eq. 2 and shown in Fig. 6(a) for run no. 1, resembles — 
the distribution measured by Jobson and deviates from the parabolic « distribution 
in the upper part of the flow. However, v, does not become constant as required _ 
for the parabolic velocity distribution, 
v, distributions for the ice- flows shown in n Figs. 6(a), and 


in the case. when logarithmic and linear s shear are assumed for a 
1 top and a bottom layer. In fact, Fig. 6(c) shows that when the top and bottom 
4 roughness are equal, v,is a maximum at mid-depthe 

_ A comparison between v, values for ice-covered flows and free-surface flows 
indicates that the values of v, for ice-covered flows are smaller than those 

for the free-surface flows in all three runs. The difference increases as one a 
moves away from the bed. The distributions of v, for the ice-covered flows _ 

are affected by the relative roughness of the boundaries. In runs | and 2, y 


distributions are skewed towards the rougher lower whereas 


on 


20 
4 
(c) 
the oussin whereas a lovarithmic velocity distribution 


symmetrical. | 


in the central regions of flow the values of v, were. constant. The present 


calculations show that this is ; only approximately true even when the r roughness 
values are same for channel bed and ice-cover surface. In cases where the 
roughness values are different, v, values are - constant even in the central 


7: 
2 
the calculated of u the stress. distribution 


a 


J=pv (16) 


‘The derivative (au/ay) was evaluated using a ‘central difference approximatio 


FIG. 7.—Shear Stress Distributions «FIG. 8.—Concentration Dis Distributions D> 
_ From Fig. 7 it can be seen that the shear-s stress distribution is very close to 
linear for all the runs. Very close to the solid boundaries, the distributions 
_ seem to deviate from the linear profile; this may be due to the approximations oy 
_ in calculating the derivative of u by Eq. 17 and the k-e model. sTtee low ob 
_ ‘After calculating the fully developed distributions ae u and v Lis , calculations 
were performed to obtain concentration distributions of a neutrally buoyant 
‘tracer by solving Eq. 6 with O = 1. After the flow has been fully developed, — 
of finite height (50 mm) ) and of uniform concentrations. “Three different vertical 
positions for injections were considered: (1) Surface injection; (2) middle i injection; — 


‘ and (3) bottom injection. The injections were either neutral, (the velocity of 


‘The 
Vimy 


the tracer was the same the ambient flow velocities at the injection point) 
a jet-type (the velocity of the tracer was twice that of the ambient flow velocity — 
at the injection point). Neutral | injections were considered for all three runs 
_whereas the jet-type injection was considered only for run no. 1. hve-sot! ie 
_ Figs. 8, 9, and 10 depict the concentration distributions resulting from three 
different source positions for run no. | with a neutral injection. The distributions 
- for run no. 2 and run no. 3 were similar to the ones shown in Figs. . 8-10. 


i In these figures the concentrations are normalized with respect to the concentra- © 


FIG. 9.—Concentration Distributions: FIG. 10.—Concentration Distributions 


for Run we. Injection) pers for Run in No. 1 (Bottom Injection) 


AG. —Velocity for Jet- «FIG. 12 —Concentration _ Distributions 


free- ited flows. The mga line of the middle injection was at a depth of 
ll cmsinbothtypesof flow. 

_ These figures show that the watuetiios in the maximum concentration values ; 
is more rapid in the free-surface flows” than in the ice-covered flows. This _ 
_ result indicates that an ice cover tends to reduce the mixing rate. The difference 

in the mixing rate is the largest when the injection is at the surface. As the | 
position of injection was moved towards the channel bed, the difference in 

Bie Figs. u and 12 show the velocity ¢ and concentration distributions for a jet- type 


a INjecuic Ol v Ss Olle emphasize 
i 
( 
q injecuon. Fig. SHOWS tne Velo aL Changes along | 
_- is of the channel from the injection location and becomes invariant with x at 7 
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ann 300 cms from the injection point. The concentration distributions charge 


in Fig. 12 indicate that the mixing is much faster in the jet-type discharge 
than in the neutral discharge and that the difference between the mixing = 

in the free-surface flow and the ice- “covered flow was not significant a 


In this paper, a procedure was been outlined to calculate the uniform, 
two-dimensional flows in channels with and without ice covers using the k-e _ 


bottom were ; made to be the same for both open-water flow and ice-covered =¥ 
conditions, thereby producing a pair of ‘equivalent flows.”’ The channel- bed 
roughness and the ice-cover roughness values were changed from runtorun. | 
The computed velocity and ¢ eddy viscosity distributions do not follow the — 

- conventional logarithmic and parabolic distributions for the whole depth of flow. 7" 7 
The velocity distribution deviates slightly from the logarithmic profile for the 
top 25% of the flow, while the eddy viscosity distribution deviates from the ' 

fam distribution through the top half of the flow. These results tend to 

{ agree with the results from recent investigations of the velocity distribution — 


_ (11) and measurements of the eddy: viscosity distribution (2), and give confidence 
to the calculations from the k-e model. 
In all cases, the equivalent ice-covered flows have larger flow depths and 
_ smaller bed shears then the free-surface flows. The resulting eddy viscosity 
_ is smaller than the free-surface flow values. However, the eddy viscosity does - 
not go to zero, and the arbitrary modification of the eddy viscosity profile 
which other flow models employ is not required in this model. 
Concentration distributions resulting from the introduction of a neutrally 
~ buoyant tracer were computed for all the runs. These distributions indicate 
reduced mixing rates in ice-covered flows compared to open-water flows. The 
difference in mixing rates was larger when the location of the tracer injection = 
was closer to the upper boundary. However, the results with a jet-type injection _ 
gave approximately the same concentrations for the ice-covered and free-surface __ 
flows, and indicate that the effects of jet mixing outweigh the difference in 
the eddy viscosity between the flows. 
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The foll b ed in thi eter ‘eens 
follo owing sym! sym ois ore in this a: 


Pes acceleration due to gravity; in Model 

= rate of of turbulent energy; 


v velocity component in vertical direction; 


longitudinal distance where flow attains fully developed state 7 

cartesian c co-ort -ordinate in longitudinal direction; 
Teg 

= mass transfer coefficient in y direction; 

Ax distance between grid points in x direction; 

= "distance between grid points in y direction; _ 
= rate of of dissipation of of turbulent energy; 


> 
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F 
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« = Von- Karman constant; ree 

= fluid density; 

empirical 

shear stress; and 
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COMPARISON OF FREE 
Free surface vortices continue to be of interest to hydraulic | engineers cers due : 
to the combination of incomplete mathematical or even physical understanding 
= of vortex characteristics with new situations in which the occurrence of vortices 
- are of concern. The detrimental effects of strong free surface vortices such 
as flow reductions, vibrations, and loss in efficiencies of pumps and turbines, 
and of structural damage, prompt engineers to design intake structures to avoid = 
“vortices. Intakes of consideration include those for ‘pumps, navigation locks, i 
7 ~ hydroelectric power, flow and level regulation, and, more recently, withdrawal 
- from sumps in Emergency Core Cooling Systems (ECCS) of nuclear power 
- stations. A fundamental problem is that it is not yet possible to analytically 
: _- Predict whether or not a free surface vortex will occur a at a particular “‘intake”’ 


for a given flow and water level. Therefore experiments, primarily — 
scale hydraulic models, continue to be used to investigate the susceptibility 
of an intake to free surface vortices. However, considerable confusion and - 
ieee exist as to the importance of scale effects in _ model studies 
due to the apparently conflicting results of past investigations. 
‘This paper reviews some of the past research on vortices and scale effects, 
and attempts to provide some order and relationship between the va ‘various studies 
as a background to a new summary of model versus prototype observations 
regarding vortex intensity. Such a comparison itself provides a critical test of 
_ scale effects due to the large change in scale and in the associated fluid mechanic © 
parameters known toinfluence vortices, 
_ Conclusions from this review and summary | are intended to help define the 7 
current state of knowledge on modeling | free surface vortices. 


Various investigators have proposed gu guidelines “for designing \ free 


he 7 , Alden Research Lab., and Assoc. Prof., Civ. Engrg., Worcester Polytechnic 
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in terms of the intake diameter 3), flow (17), velocity (6, 12), ora , withdrawal ~s 
Froude number (26,25) without explicitly including a measure of circulation. _ 
Basic guidelines, such as by Anwar (1), require explicit knowledge of the = 
circulation, a factor not known a priori. The acceptable versus unacceptable 
pumped storage ‘intakes presented by Pennino and Hecker (22) were tested — 
against the available practical guidelines with essentially no correlation, particu- 
i - with regard to the concept of a required submergence. The inefficacy 
of the available practical guidelines has also been mentioned by Dhillon (8), 
§ 4 and such guidelines are not intended to apply in cases of strong local circulation * = 
~ as generated by offset intakes or screen blockage in ECCS sumps. Cas _ 
- Despite decades of work, therefore, the present state of knowledge is such — + i 
that reduced scale hydraulic model studies or some other | experiment is required 
_ for each case where it is important that no strong vortices occur. A possible 
ss with scale models, however, is that free surface vortices are subject 
to prediction errors due to the impossibility of reducing all pertinent forces 
_ by the same factor. Although the predominating inertial and gravitational forces 
are reduced similarly in the Froude-scaled models used in such studies, viscous Ea 
and surface tension forces cannot simultaneously be reduced as much, and { 
the extra influence of these forces on modeling air-core vortices is called a 
“scale effect.’’ Considerable research has been conducted on the magnitude 
of scale effects in modeling air ¢ core vi vortices and on developing a 


@ « with other flow phenomena subject to viscous scale effects, it is important — 


that the model Reynolds number be sufficiently high. Both Anwar (1) and Daggett 
and Keulegan (4) give a minimum suggested Reynolds number t of about 3 x 
—10*, based on inlet flow and submergence or intake diameter, respectively. 
These investigators indicate that viscous scale effects on air core vortices are 
absent when the Reynolds number is above this limit, thereby implying that 
_ models satisfying this criteria can be operated at Froude-scaled flows . That — 
oe is a sufficient criteria has not been universally accepted, although an early 
study by Quick (24) s suggests Froude scaling is s adequate as long as the Reynolds ig 
number influence on vorticity generation is taken into account. A recent a 
by Dhillon (8) indicates good model-prototype agreement of vortex strength | 
for a 1:20 scale model, satisfying the Reynolds number criteria, operated at 


Testing at higher than Froude-scaled intake velocities (flows) to overcome a 

Z ‘scale effects of modeling air core vortices, including operating models up to ial 
7 prototype intake velocities, has been proposed by various investigators. al 

such tests maintain the geometrically scaled water depth, such that approach | 
flow velocities are increased in the same proportion as the intake velocities. 

7 Haindl (13) and Jain et al. (19) reached such conclusions based on experiments a 
_ using various sized circular tanks, while Denny (5), Iversen (18), Dicmas (9, 1, 

- Zajdlik (28), and Chang (3) have reached the same conclusion using various _ : 
_ sizes of a particular pump intake geometry. Both types of studies show that 
_ the required increase in flow above Froude scaling is related to geometry and 
Telative | submergence, such that a constant factor for - flow increases ead 


q 
| 
| 
| 


a Hecker a 1). . Conclusions derived from such studies regarding the usefulness 
of testing at higher than Froude-scaled flow are based on a limited change — 
in facility size, typically by a factor of about 3, although larger changes in 
scale, such as by a factor of 8, were achieved in some cases using rather 
small outlet diameters. Denny and Young (6), based on one observation of | 
a; prototype ‘pump previously modeled at a 1:16 scale ratio, also’ concluded 
= using prototype velocities in models of such a scale ratio may be proper. 
Fo some scale ratios, the use of the equal intake velocity concept would seriously 
Barer the primary Froude-scaling criterion used to achieve proper appre 
~ flow patterns and the resulting circulation at the intake. A recent series of 
ry model studies by Dhillon (8) indicates that use of prototype velocities” 
in “models produced highly exaggerated vortices, incompatible with prototype 
observations. Similar were apparent in an earlier study by Linford (20). 


predicted by the model to a relatively. higher prototype value, has the inherent 
“difficulty that ambient circulation may change with geometric differences at 
higher elevations, and such a technique would, therefore, be most applicable — 7 
to two-dimensional geometries. In addition, other studies have shown that the ‘ c 
required correction is a function of the basic geometry involved and of the 


a) and Quick (25), is often not evaluated, whesens this mechanism is heen evo 
cases in which no boundary discontinuities exist. 
 Aclear that the vortex air core characteristics are not affected 


- various degrees of surface tension but constant viscosity. Daggett and Keulegan = 


_ indicated that flow characteristics (discharge coefficient) varied little for a 
threefold variation in ‘Surface tension, although the relationship between the 


by Hughes (15), ‘who also provided | some rationale for the constant velocity 
modeling concept by considering the effects of surface tension on air core 4 
vortices. Jain et al., showed that the critical submergence for incipient air — 
- entrainment was essentially constant when the surface tension was reduced — 
by about 50% for any given value of intake velocity, outlet diameter, and ¥ 
_ However, the range in Weber number gi given in both references comes mainly 
from changing the velocity and outlet diameter, either of which also changes 
iy the Reynolds number and, therefore, the vortex core. Isolation of Weber number 
effects would d require plotting a vortex-related parameter versus the Weber ee 
for constant values of circulation, geometry, and Reynolds number. Such a 
plot, however, would be difficult to generate due ae the stated interrelationship . 


i not been found. A generalized approach to such flow increases and a technique = ; 
q 
— 
: considerably higher than dictated by Froude scaling would tend to dissipate _ ; 
e conc onose Tai al (19) of correcting the critical cuhmergence 
and the VOrex all COre Charactcrisucs Was 


 /ahe as compared, for example, with the maximum tangential velocity and = 
: the radius of the air core tip, respectively. Since these latter parameters may 
= not scale directly, the difference between model and prototype Weber numbers 
_ may be difficult to determine. The recent analysis of Yildirim and Jain Qn 
indicates that the effect of surface tension on air core vortices may not be 
negligible, and, in view of the above comments, the question of surface — 
Irrespective of the relative influence of surface tension and viscous hh 
= the air core, it seems evident that vortices with air cores extending near 


: = or into the inlet are more subject to scale effects than surface a 


And corresponding subsurface swirl which can be made visible by injecting 
dye into the dimple. In a full-size installation, due to the absence of scale 
effects, the transition from a dimple and coherent subsurface swirl to an a 
core will be relatively more rapid with increasing swirl intensity compared with _ 
a model. And the initial condition of coherent swirl from the surface to the 
: a dye core vortex, would be a reasonably safe limiting condition q 
cases in which air core vortices would be detrimental. Based on experience 
| the Alden Research Laboratory, the delineation of dye core vortices as part 
the normal model testing program is notadifficulttask, 
Given the small change in size used in most experiments to define the required | 
_ flow increased above Froude-scaled flows, and considering the general neglect 
of such factors as the transient nature of vortices, the relative effect of boundary 
- layer vorticity compared with other sources of vorticity, and the effects of - 
ambient flow turbulence on vorticity, the concept of increased model — 

— above Froude scaling could be further examined. However, the past emphasis 
- on quantifying scale effects of modeling air-core vortices is somewhat perplexing 

in that few final designs developed by using scale models have such strong 
vortices. Perhaps a more pertinent concern is the reliability of models as predictors 
of weak vortex activity. For example, does a model dye core filament and 
surface dimple at Froude-scaled flow represent an air core vortex in the prototype? 
This matter is addressed by comparing information on vortex activity from 
model studies with similar for the prototype. 
‘“ The ultimate test of the ability of a reduced scale hydraulic model to predict d 
free surface vortex characteristics is to compare model data to corresponding 
prototype observations. Because of the inherent difficulties with field observations 7 
and the influence of other factors besides scale effects which can . contribute: 


to differences between model and prototype vortices, a significant number of 
_ such comparisons are needed before general conclusions are warranted. 
_ Included within the first category of inherent difficulties with field observations 
_ is a natural tendency for an observer to be impressed with the size and intensity _ 
of prototype vortices, and a belief may develop that they are” more severe és 
than those predicted by a hydraulic model. This belief may be partly due to — 
‘the increased levels of ‘surface turbulence and noise, which undoubtedly do 
not scale directly. The difference in time scale also makes it appear that the 


Deena vortex is more persistent, whereas the prototype may, in ‘reality, — 


| 


_ is also difficult to » diatinguich mec various strengths of vortices in the field, i. 
and the depth of the air core cannot usually be determined. Other factors wae — 


1 attention to and boundary roughness 
“ the model, aged approaching the inlet has a first order ef! fect on the resulting 
eal 2. Viscous bigs ef fects on modeling appurtenant flow devices such as screens, a 
baffles, vortex suppressors. 
“Minor” prototype topographic or structural chenge. from that tested in 


cofferdam cells, and wing walls. = 
5. Ambient ent density stratification. palit 
number of model- -prototype of intensity have been 
= published, and, as sufficient information was made available for an independent — 7 
reevaluation, the following will be included in this summary. A study by Linford — 7 
= (20) using a 1:200 scale model of the Kariba Dam and powerhouse intakes ~ 


» 


=F showed that air core vortices were produced in the model when operated near 
a 
a 


prototype velocities, whereas no air core vortices were observed in the fi field. 
Reasonable c comparison of vortex - strength was achieved at near Froude- scaled 
flows. Nelson and Johnson (21) published photographs comparing model and 
prototype vortices at the upstream intakes of the Eisenhower and Snell navigation | 
locks, and concluded that the prototype vortices were larger than in the model. 
a conclusion i is not, however, obvious from the photographs since the 


> A similar « comparison for the Demopolis lock, by Nelson and Johnson (21) n 
and The U.S. Army Corps of Engineers (16) appeared to give a reasonable _ 
correlation between model and prototype, and this was confirmed by viewing — 

- movie with comparable scenes of both the model and prototype vortices. 

i lock models were all operated at Froude-scaled flows. Recently, Dexter 
and Ziegler (7), using a 1:120 scale model of the Grand Coulee Third Rewer 

_ Plant, predicted weak vortex action at the penstock intakes whereas no ee 
have been observed in the field. This comparison is based on model flows 
of about four and one half times the Froude-scaled flows which occurred during 
the prototype observations. However, there are questions as to whether the 
= flow circulation was properly simulated and about the detail of the 
sg observations, | so the apparent exaggeration of the model vortices has not 

To obtain the e larger. data base required for a reasonable analysis, questionnaires. 
requesting specific information on projects for which model-prototype vortex 
Bee | was observed were sent to 65 organizations, both domestic and foreign, 
_ which have had d occasion to be concerned with vortices from a design, testing, 
or operational perspective. ow! | hoqmeg 
‘The returned questionnaires were reviewed to enduin’ which provide useful © 


_ model-prototype comparisons of vortex activity and information on prototype 


operating or = known thereof, while vortices occurred. The 


| 
/ 

| | 

q | 
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‘Bear | Creek flood control 


outlet (TVA) 
Bear Swamp 


storage, upper structure 


England Power 


upper 
Service al 
Colorado) abd 
Dardanelle Lock (USC 


mA, 


4 Demopolis 


of EY 


Dover Dam of E) 

~ 

Eisenhower Lock (usc 


of E) 


Ffestiniog pumped storage; 


lower reservoir structure 
-Foyers pumped 
upper reservoir structure 
(Scotland Hydro Board) 


Grand Coulee Third Power 


| 
| 
Heart B Butte Dam (USBR) 


Jocassee pumped storage, 


structure (Duke 


pines 


single sluice gate, 

cylindricaltower 

single 40 ft diameter open 
intake; V = 9 ft/s, S =a 


upper; multiport side imake 
~V~ 25 ft/s, 
oat vee a 
multiport (8) at upper gate 


18 adjacent sluices, 

O<S<S4f, 

0<V<47ft/ 
wpe 
four horizontal intakes, 32 _ 


aa x 22 ft, S$ = 40 ft, V= 


four openings, 
total 68 ft x 24 ft, 
«23 ft, 4S (4 


Six adjacent penstock 
_ intakes, each 33 ft x 45 


slow circulation, F, 


TABLE 1 Comparison 


(3) 


= 1/20F, =l 


intermittent, some pulling. 


air bubbles at F,= 1, L, © 
PIG | 
=i, 


=1, 


swirl, F, = 1, 

i 
no vortices for F, = 1 0. 


V,) L, = 1/60 


depression at 2 < <4, 


dye core no air 
entrainment, 1/120, 


V = 20 ft/s, S = 180 


morning glory spillway, v i 
3 ft/s, 0<S< 54ft 


two tower intakes, V = 6.5 


negligible, L, 


small vortices between gate : 
guide piers, L, = = 1/50, OF ie 
F=1 


drawing vortices for V. 


a 


L, = 1/33, 
ingen 

= 

| 


of Vortex Intensities =: 


+ 
4) 
_ depression. 


r intermittent, sometimes 3 ft 


ameter, pulling air and 


no almost no motion 
tnt 
two audible vortices, up to 
8 ft diameter intermittent 
tocontinuous 
none at normal 


_ submergence 


turbulence, also 
hazardous 


3890 


no problems a 


swirl depression 


‘no 
ior i: ned 


_ Surface turbulence, no ne 


small, | in.-2 ee. transient 


vo ices 


no no probleme. 


Operation 


no problems 


air entrainment, possibly 
increased lock chamber 


— Comments and Key® 


pos correlation of type 
and size (b) 
prototype vortices seem 
more frequent and stable, : 
model retested 


(comparison based on 
= 1); see Fig. (a) 


concern for ice entrainment 


led to conservative design 
— 
vortex free 
_ developed using model 
model study after prototype 
_ problems, checked initial 


design (16,21) (b) 
en ( ) (b) 


larger (21) (a) 
near equal velocity — 
exaggerated vortex 


4 


model (conditions at F, = 


| 1 used in Table 2) (b) | 


good comparison at model aa 

velocities somewhat more 

*, than Froude scale (b) 

comperisoa based on Units 
19, 20, and 21; more 
severe model vortices 

maybe caused by 


simulation (c) 


| comparison based on final 


design with six radial we 
crest piers to minimize 


vortices 


SURFACE VORTICES 1249 | 
- 
| 
| 
a 
problems" 
| 
“surface dimple at S, = 20 ft | 
| 


19810 


Kariba — pie | six intakes in borderline vortex at F, = 
scheme (Zimbabwe slopingright air core at > 3, L, = 
Zambia) embankment, each 54 ft 1/120 =" 
x 32 ft, V = 3 ft/s, 50 ft 
Ludington pumped storage, | six adequate horizontal nimal wh with dim, dimples 


upper structure inlets; each ~30 ft x 35 
(Consumers Power). ft, V = 11 ft/s, od 
Muddy Run pumped four tower = vortices (no air) above 
Storage and upper 470; air 
structures (Philadelphia = entrainment at elevation 
“a | 466 or less L, = 1/35 
ae 
single horizontal inlet at _ negligible, L, = 1/46, F, 
‘upper structure end of channel, V = 5 
(Northeast Utilities) “ft/s, S = 69 
Oroville Dam diversion two adjacent 35-ft diameter strong, ‘noisy, air core, 
tunnels (California 73 tunnels near reservoir | ¥ sucked scaled trees, __ 


Department of Water at bottom, 0 < S < 690 ft, 0 f persistent at S < 100 ft; 


Resources) ds transient at S > 100 ft, 
ier 18 ft 14 air entraining, stable as= 
auxiliary intake (India) $< 47 ft, V= 11.5 23 ft, L, = 1/40 
Reacteurs G2-G3 Marcoule a openings in sloping _| large noisy acrated vortex ; 


cooling water intake embankment, each = 5 
Ca 
Snell Lock (USC upper gate 


of JSC of E) 


a ‘Taum Sauk pumped shaft, — to D, air entraining vortex at 
storage, upper reservoir | = 46 V= = elevation 1,508 ft; less 
intake (Union ion Electric severe at 4% 


i 

Treasure Island Pumping two vertical suction air = 


ation (USC of columns, V = 2.4 ft/s, S| 1/12,F=1 | 


Note: 1 ft = 0.305 m; | in. = 25.4 mm; 1 ft? = 0.0283 m’. 


| TABLE — 
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HY1IO 
Continued 


swirls (initially had 3-ft 
_ diameter vortices during 
Start-up) ice 


4 vortices as level approaches 


elevation 470 ft 


air core for 100 ft< S< 
230 ft, passed trees, 
strong audible vortices at | 
uncontrolled flow 


large noisy aerated vortex 


vortex with depressed cor core 

air entraining vortex at 

a elevation 1,515 ft or less” 
air entraining vortex, core 
diameter = j 


FREE SURFACE 


near equal velocities 
7 exaggerated \ vortex in 

model (comparison near 
m4 F, = =1 with Units 5 and 6 
reservoir relevation 
1,549 ft (20) used in 
perforated wall used to 


dissipate two strong 


nad? 


objectionable turbine noise 


at extreme low water, <_ 


| elevation 470 (operate at 


blocked racks dueto 
suction of debris 


boa 


run turbine at 1 /2 flow ow at 

low surface elevation to 
avoid rough operation — 


| vibration and adverse a 

_ effects on machinery 


model vortices; initial 

vortices in prototype not 

observed since start- ee 

2 

model predictions based on 


prototype 


persistent at higher 
submergences; otherwise 


agreement 


model tests after prototype 
problems; remedial 


: developed using 


prototype results (b) 
prototype vortex seemed 
larger (21) (a) 


plan was to use floating 
grid in prototype, if 


model to prototype my 
pari: is son Froude 


or or (floating grid) 


evil 


no problems (including 
initial start-up) 
| 
n 
&§ 
‘ 
. | — 
| 
develop vortex 


= Unfortunately, a a large portion of oe returned questionnaires \ were not useful 
_ for a variety of reasons, including provision of information on a model study — 
only, use of imprecise or contradictory terms, or generally insufficient informa- 7 
tion. Follow-up inquiries were made regarding projects of interest and other 7 
Projects previously in various publications. The numerous reports, 


The information on some projects is considerably 1 more » detailed, whereas the 

7 filled-in questionnaire was relied on for other cases. It is consequently important _ 
3 this review to concentrate on basic trends rather than detailed observations —_ 
 Companteons in which the final model design and the corresponding prototype 
had only weak swirls and no vortices were considered meaningful in that they © 

showed field conditions were not worse than expected. 

‘ Projects for which an adequate comparison between model and field vortices 


aes are listed in Table 1 which also gives a basic description of the intakes, 


BLE 2.—Number of Projects with Indicated Model- ‘Prototype Vortex Comparison 


(a) Prototype vortices stronger or more 


1 
q 


‘(b) Model and prototype vortices 
(c) Prototype vortices weaker or less 


— 


information ‘regarding: the model tests, and any additional pertinent 
_ Whenever possible, the wording is that used on the material supplied, although 2 
some changes or interpretations were made for consistency and to concisely 


intakes, e. g., ‘for navigation locks, pump intakes, or screen n wells, outlet control 
- structures and sluices, specialty intakes for diversion tunnels and water-supply- Pe 
‘regulating tanks, cooling water intakes for steam electric stations, and a consider- 
number of pumped storage intakes. 

__ A basic summary of how many projects fit the case of prototype vortices 
stronger than in the model (a); essentially ‘comparable vortices (b); or vortices 
less intense in the prototype than in the model (c), is given in Table 2. Since 
a comparison is most when related to any model testing techniques 

i to overcome scale effects, this summary is provided according to the scales 

: flow rate used in the model tests. The major observation apparent from Table 
2i is that most modeled vortices compared reasonably well with observed prototype 
Five projects seem to have stronger or more persistent Prototype 


: S information on operating problems was considered important since it is relevant nyt 
a 
sib: anal 
“See Ffestiniog and Kariba, Table 1, and Dhillon(9), 


had Reynolds numbers equal to or greater than the seneaianiial minimum 
values given by Anwar (1) and Daggett and Keulegan (4). Since these projects 
shad occasional air " core vortices, some scale effect for Froude models of air — 
core vortices may / be indicated. That this is a weak conclusion is indicated | 
A by the following comments. Two of the five projects are the previously mentioned 
Eisenhower and Snell lock intakes (21) for which no conclusive comparison © 
could be made. A third case refers to the apparently more persistent prototype 
vortices at the Oroville Dam diversion tunnels at high submergences during 
a flood. However, few model observations \ were made | on sa once activity at 


= the diversion “phase, and a ‘comparison of persistence is not based on — 


A STRENGTH LEGEND 
‘SWIRL = | 
F 
DIMPLE 
DIMPLE WITH DYE cone 


PULLS BUOYANT DEURIS 
AIR CORE TO INTAKE b) 


For Prototype, Determined by 
For Prototype, Determined by 
Qhange in Nowe Type and intensity 


ATEX STRENGTH 


‘OR 


100% = FROM 
EL 1599 TO 1560 ——» 


— con Froude Sealing) 


At the Ludington upper reservoir not predicted 
te model occurred only during initial operation of two out of six units. No — 
coherent vortices have been observed in the field for years subsequent t to _ 
_ _ This case raises an interesting point regarding scale selection and the minimum _ 
Reynolds number required to avoid scale effects. Without the perforated wall 
of the final design, two strong and persistent air drawing vortices existed in 
_ the model. Suppression of these vortices with the perforated wall, or with any 
other similar technique, involves dissipating angular momentum by form loss, 
as well as reorienting flow patterns. The former factor may be dependent on 
the Reynolds number based on the velocity and size associated with the energy 
5 - dissipating element. It is possible that excess energy and vortex dissipation 
~ occurs in small-scale models in which vortex suppressors are used. When selecting | 
a model scale ratio and when testing vortex suppressors, this should be . 
consideration as should the withdrawal Reynolds 


The last case hin category (a) is the Bear Swamp | ‘Teservoir 


q 
CW tae PERCENT TIME INDICATED VORTEX OR STRONGER EXISTED 
q 


Project and Type 
Beach City Dam 
flood control outlet 
 (USCofE) 


Centrale 
hydroelectrique de 


Nationale du 
Rhone) 
Centrale Thermique 
de Cordemais, 

pump intake (EDF) 


Enid Dam flood 
control outlet q 
4 (USC of 


Glen Canyon Dam 


198: 


OCTOBER 1981 


* Intake 
~ 


towerinopen 
reservoir, six 
adjacent gated 
openings, each 7 ft 
x 0< Sq< 
40 ft, V varies 
withhead 

single bulkhead shaft 
for tens tunnel” 


intake to six 
turbines, v =5 
ft/s, S = 26 ft 


vertical pump 
column, 16 ft 
diameter, V = 3.3 
ft/s, S= 8.2 


tower in open oa 
reservoir, two 8-ft 
16-ft intakes, V 
= 30 ft/s, S = 43 
open approach to 24 
bellmouth, V = 9 
‘ft/s, 3h<S< 


_ Kori Nuclear Station 
Cooling Water 
Pump (Korea 
Electric Company) 
Mohicanville Dam 
flood control outlet 

(USC ofE) 


four vertical pump 


columns, bell 
_ diameter = 7 S ft, 
V=82ft/s,S= 
tower in short 
 forebay, three 7 ft 
x 12 ft gates, no 


conduit, 0 < Sg < 


_ Nimbus Dam power 
intake 


Fal 
k (USC of | of 


25 ft, V varies 
| head 
open approach to 41 
powerplant intake 
rm in dam face, V = 

3.6 ft/s, Sg = 10 
gate sill, V = 6 
ft/s maximum, S 


Prototype vortex 


ft-3 ft diameter 
vortex at 2,200 


no direct 


observation, 
inferred from 
model study 
no direct 
observation, model 
study showed hil 
intermittent air a 
drawing vortices 
nodirect 
observation, model 


| 


study showed 
intermittent air 
drawing vortices at 


minimum 
3-ft diameter 


full gate opening — 


Operation 


Operation Information is 


Comments 


Set 


for largest vortex, | occurs at particular 


wake 


dislocation and 
grating 


irregularities = large probability that 
approach flow problems are due 
resulting in to prototype 
reduced vortices” 

large probability that 


considerable noise 
aa problems are due 


and vibration 

(judged dangerous to prototype 

to support vortices 

Structure) at 

minimum 
submergence, 

depriming and unit 
trip 
reduced efficiency, 
possible gate 


flow affects stilling 


swirling | 


submergence, two 
transient vortices: 
within intake, | ft 
diameter, 6 ft core 
depth 
nodirect 
observation, 
inferred from 


study 


mini eddy near right 


5-ft diameter ha 


vortex 


ing! 
— 


audible 


basin performance 


No problems 


large probability 
problems due to | 


prototype vortices Ay 
occurs at particular a) 


gate opening and 
submergence 


rough operation of raft used to dissipate > 

Operation achieved — 


th Bar 
subsequent to field sit 
resulted problems, model 


in fatality ‘Study was 


initiated; 
comparison with 


excessive noise and 
vibration, 
contributed 
gearbox failure 
| 
slight vibration 


hazardous flow 


TABLE 3.—Additional Projects with Vortices for Whic pa 
ty 
_ 
| 
= 
— 


 Shadehill Dam Morning Glory 2-ft diameter, audible | no apparent ill similar to Heart 
(USBR Spillway V~3 | vortex; pulledin | effects on Butte Dam 


ft/s,0< S< 40ft ice sheets and 50 structures or Spillway, but 


4 


= 


operation, no 


structure 


incorporated in 


intake, for which considerable model and prototype data were recently made © 
z _ available by Castro (2), including the percent of time (persistence) that vortices : 

of various strengths occurred. The method of differentiating between vortices. 
of varying strengths, which was s developed at the Alden Research “prec omagll , 
a is indicated in Fig. 1. The vertical penstock of this project has a single round 
opening at the end of a short approach channel near the edge of the reservoir. — 

Basic results of carefully retesting the model, in conjunction with prototype 
observations by the same person, indicated that prototype vortices were stronger = 
and more persistent than occurred in the model at Froude-scaled flows. Fig. S 
1 shows that flows of about two or two and one half times the Froude-scaled - 
flows were required to reproduce the field persistence observations. _ oro 

Not all aspects of vortex characteristics were reproduced, however, with 

the increased model flow. For example, the trend of vortex intensity versus 

reservoir elevation was not identical despite modeling the entire upper reservoir, 
‘indicating that factors such as relative differences in friction for the > reservoir a 

flow may have changed the generation of vorticity. 

Of the 16 projects in category (b) showing good model-prototype comparison, — 

two compared well when the model flows were increased from 2-4 times Froude 

scale, and the remaining 14 comparisons all being based on Froude-scaled flows. — 

_ Of the latter group, 11 concern good comparison of negligible or weak vortices, 

a significant finding in that negligible scale effect for weak swirls or dimples 
S implied. The remaining three cases concern good comparisons of air a 


vortices, and these cases are not as instructive in that any scale effects would 
dictate that an air core model vortex certainly represents an air core prototype 
a ‘The single case in category (c) is the previously discussed tentative conclusion © 
regarding the Grand Coulee Third Power Plant (7). 
“7 In general, final intake designs which were developed from Froude-scale model _ 
tests to be vortex-free were indeed vortex-free in the prototype, and ad is 
that had weak vortices in the model had weak prototype vortices. In particular, a 
there is no reported case in which a negligible model vortex corresponded a e 
a sufficiently strong prototype vortex that it produced operating problems. The _ 
operating problems to be discussed below were either predicted by th the m model | 
or were for projects for which no initial model study was conducted. — ee Sniviegs 
In addition to the in Table 1 for which information on 


| 


J 


on projects which lacked model- -prototype comparison of vortices. These various S, 
7 additional projects and a brief description of the intake types and the operating 
problems, or lack thereof, are listed in Table 3. A summary of project operating wa 
| conditions related to vortices, from both Tables 1 and 3, is provided in Table 4. 7 
It does seem ‘that ‘prototype vortices have produced sufficient “operating 
“ _ problems that reasonable measures should be taken to avoid strong vortices. — 
Of the 20 projects which are known to have (or have had) reasonable sized 
_Vortices, 14 have experienced the various types of problems listed in Table | 


4. A considerable amount of operating dif difficulty i is associated with prime movers, A 


particularly for “wet-] -pit”” pumps. ‘Pumped storage projects for which air core 


vortices were observed in the model at low water levels produced unacceptable 
; ‘turbine operation in the field at corresponding conditions, thus requiring flow | 
reductions to avoid the vortex condition (e.g., Taum Sauk and Muddy Run). ; 
‘Air entrainment may also be a problem in navigation locks in which the release — 


; ~~ 4.—Summary of Projects with Information Regarding Effect of Vortices on 


AP AL Prototype Operating Problems* Number bus 


Vortex reduced flow, or flow reduction to avoid vortex 
Vibration orpump) deol 
; _ Other pump or turbine operating problems (i.e., rough opera- | seend 
tion, objectionable noise, lower efficiency) pe, 2 
Detrimental effects on stilling basin performance _ 


air in the chamber excess ‘turbulence. Also” to be noted is 
six projects have experienced no problems while vortices occurred. 

The p possible occurrence « of vortices has generally defied analytic prediction, — 
and experiments, usually scale model tests, are needed to ensure vortex-free | 
designs. Such scale models should have Weber and Reynolds numbers above | 
published critical values to minimize scale effects(1,4). 
To help | clarify apparently contradictory information previously published 
regarding scale effects of modeling free surface vortices, comparisons. were 

_ made between model and prototype vortex characteristics of various intake 
_structures. Since scale effects axiomatically have only a secondary influence, at 


a large change in scale was considered desirable to make any such effects 


| 
= 
i 
| 
6 
} FOr each project, the major problem, if any, 1s selected and no project is listed more a ‘ 


a flow pattern. The concept of equal velocity modeling seems relevant only to — 


‘more obvious. The cannot encompass all available 
information, but together with other previously | published data, provide lag 7 
to reach the following 
4: Sufficient operating problems due to vortices have occurred to warrant — 
pent efforts to minimize vortices at intakes of various types. _ ates 
a 2. Hydraulic models operated at Froude-scaled flows to predict vortex intensity 
s some scale effects when ‘simulating air "core vortices. Compensation for 


- scaled values, but large it increases are e not justified and may distort the approach ; 

me large models for which small flow increases yield prototype inlet velocities. _ Sil 
3. Negligible scale effect was evident in models with Froude-scaled fiows 
= predicted that only sw swirls and surface dimples, but no air core vortices, 


would occur in the prototype. A possible exception involved a vortex dissipator 


‘itself subject. to scale effects. ne Wedel 


addition, the following recommendations are made: 
1. Considerable attention must be given to scaling the approach topography B 
boundary roughness so that the “vorticity flux approaching the intake 
cosrectly modeled. The difference in Reynolds number and friction factors 
_—— model and prototype boundary flow should, therefore, be evaluated. 
2. Model tests should include i. careful delineation of various vortex types, 
¥ including dimples and dye cores, and notation shc should be made abet type, 


location, and, in in particular, persistence. 


- acceptable and unacceptable vortex conditions since this represents the initiation 
7 ¢ of coherent subsurface swirl. Since modeling such a dye core and surface dimple 
M _ involves essentially no scale effect, models can be operated at Froude flows — 
to detect such | flow patterns. Dye cores should be present for less than 50% 
4. When vortex suppression devices are used in scale models, consideration 
should be given to possible excessive energy dissipation due to low model Reynolds 
numbers based on the characteristic veloeity, and dimension of the dissipator. i 
‘This effort was initiated to assist in meeting the € objectives | of the ASCE 2, 
Task Committee on Intake Vortices. The information and comments supplied r 
H. O. Anwar and G. Young of the committee are appreciated. 
Considerable assistance and information was supplied by a number of organiza- _ 
tions listed Table and help of D. of the U. S. Bureau of 
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F Froude number ratio, model to 
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TECHNICAL Notes 


as yet, to the point where they warrant publication as a Proceedings paper in a Journal, a ; 
publication of Technical Notes was authorized by the Board of Direction on October 16-18, — aie 


= . An original manuscript ay two copies are to be submitted to the Manager of 9 = 


ws Professional Publications, ASCE, 345 East 47th Street, New York, N.Y., 10017, along a 
with a request by the author that it be considered as a Technical Note. © eet eS ooo 
2. The two copies will be sent to an appropriate Technical Division or Council for review. 
 @ If the Division or Council approves the contribution for publication, it shall be returned 
to Society Headquarters with appropriate comments. = = j= | 
_ 4, The technical publications staff will prepare the material for use in ‘the earliest possible 
issue of the Journal, after proper coordination with the author. 
» 5. Each Technical Note is not to exceed 4 pages in the Journal. As an approximation, — 
each full manuscript page of text, tables, or figures is the equivalent of one-half a Journal 
6. The Technical Notes will be grouped in a special section of ‘each Journal, 
4 7. Information retrieval abstracts and key words will be unnecessary for Technical Notes. 
_ 8. The final date on which a Discussion should reach the Society is given as a footnote .- 
each Technical Note, | 
9. Technical Notes will not be included in Transactions. — re 
4 10. Technical Notes will be included in ASC E’s annual and cumulative subject | and author ‘ 
indexes. 
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2. The author’s full name, Society cmbnate grade, and a footnote reference stating present oa 
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- The letter symbols used must be defined where they first appear, in figures or text, and arranged — 

alphabetically in an Appendix.—Notation. APA a 
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published by the American ‘National Standards Institute and to the Authors’ Guide to ‘the 

6. Tables must be typed double-spaced (an original ribbon copy and two duplicate copies) _ 
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UNCERTAINTIES IN WATER DISTRIBUTION SYSTEMS = 
Derrick A. Clarke," Edward A.McBean> 


"and S. A. Al-Nassri,’ M. ASCE 
Water networks requ require at major of the investment 
of water works. As a result, over the years the analysis of water distribution 
networks has been the interest of many researchers. While these activities resulted 
in the development of general computer programs which can handle the analysis 
| a large network in a matter of seconds at nominal cost, in contrast, little — 
work has been done on the uncertainties associated with the analysis of water 
distribution networks. Hudson (6) states that 
_ Unfortunately the design of water distribution systems, based on poor _ 
_ a data, looks just as good in the computer printout as a design based on vo 
most accurate data that can be obtained. 
As a means of examining the uncertainties, the objective herein is to examine _ 
“the sensitivity of the pressure difference, AP, in the analysis of a distribution 
“network as a result of uncertainties in the parameters involved in the basic 
- equation describing the flow. _ The First Order Analysis method developed for 
_ evaluating the pressure sensitivity is used to assess the required pressure — 
of the give hes been im service, thes the 
Hazen-William’s formula can be to express AP in terms 
in the 1 model error term which i is for additional 
me differences due to fittings within the network; y = the specific weight; 4 
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Q = the flow; n = 
the — = the Hazen- William’ s and D= 


conclusions to follow are equally applicable to any other equivalent formula. Big 
The application o of first order analysis © Eq. 1 gives (see Refs. 3 or 10 for 
Suse 


+ (log, 

in a which the bar (~) in Eqs. 2 pie 3 = the mean value of the east 

Q = the coefficient of variation of the parameter; and p, = the coefficient - 
of correlation between the parameters noted as subscripts.§ = 


The last three t terms in Eq. .3 account for correlation between ‘the parameters - 


a variance of | the parameters in n the basic equation describing the flow - 
each pipe and the model error are known, Eqs. 2 and 3 can be used to evaluate > 
the mean and coefficient of variation of the pressure differential AP. 


- In order to illustrate the application of the. methodology of first order rarer : 
_ (FOA), a simple, 1l-pipe network is considered. The sensitivity of pressure “ 
7 differences, AP, for a specific pipe in the network will be examined, as a 
_ result of uncertainties in both the eT. involved and the mathematical © 


5 Consider, e.g., one f pipe of the 11-} -pipe ye network, 305 m in length and 304.8 

mm in diameter. For a pipe that has been in service for 20 yr, the mean value ~ 
7 of Cyw is 105, and the deviation, is 10, as obtained with the 


. 2 ‘After 20 yr in in service the diameters of the Pipes: are © subject t to change due 


as s it is very much a function of the quality « of ‘the water surveyed, the filtering 
efficiency in the water treatment facility, and the design of the network system. | 
_ However, from the writers’ experience, a 5% change is reasonably normal, 
. will give for the pipe under consideration a standard deviation op) of 
The variation in the ideal of the baitindiie. slope of the Hazen- Williams 


- formula (Eq. 1) is assumed, for the present problem, within the range 1.82-1.88, 
_ with a mean value of 1.85. [This assumption is on the conservative side and 


rough, of the Moody diagram. Moore (8) introduced a modified Moody diagram 
with lines: of constant exponent ‘‘n”’ varying from 1.7-2.0.] Fora given relat relative | 
roughness e/ D and range of flow, n’’ can be eeasily obtained. 


ic slope; Lf 
the internal 
been utilized 
= 
| 
J 


“all the fittings (e.g., valves, fire hydrants, and bends). The additional head 
loss can be accounted for by considering an equivalent pipe length being created — 
_ by these fittings. The recommended fittings (7) for the pipe under consideration — 
' are: 2 Gate Valves, 3 Tees for fire hydrants, and | Tee piece for branch mains. : 
x The equiv alent lengths « of pipe resulting from the presence of the above fittings = 
_ were determined from the table produced by Hammer (4), and found to be 
10.4 m, 18.3 m, and 6.1 m, respectively. Therefore, the total equivalent length — 7 
is 34.8 m and the model error, A,, for the 305 m long pipe is 339.8/305 = 
1. Assuming 10% coefficient of variation, 
Values of 0. 67, —0.64, and 0. 67 were obtained for the coefficient of correlation | 
between Cun ond | Q, Cu nose and D, Q and D , respectively. The coefficients 
of correlation between both Cyw and Q and Q and D were determined from 
_ the values of C,w, Q@ and D obtained from the variation of C,,y, for the pipe, | 
using a standard computer program for network analysis (2). The preceding 
second coefficient of correlation (—0.64) was determined d directly from tables — 
established by Williams (9) assuming an age of 20 yr. im 
‘Using Eqs. 1 and 2 and the aforementioned coefficient values, a mean value 
of 18.5 KN/m* was obtained for the pressure differential, AP, across the pipe, ; 


-= a coefficient of variation, 2,,, of 0.337. This implies a standard deviation, 
Crrrenia FOR Dererminins Reauireo “Accuracy IN PRESSURE Measurement Devices 
oer the results obtained using FOA principles, the following a are 
for the determination of the required instrument measuring 
‘Pressure i ina a field test testing program: 


oa If the standard deviation of the pressure differential, o,,, obtained from 


the length of time the pipe has been in service, then the measurement devices 


OA is less than the additional pressure differential, A(AP), resulting from 
Tf: < BAP ron. 


7 
E + Car, < AP JFOA 


which E = the expectation operator; and = the standard deviation 


the measurement « of the | pressure differential between points A and B due to 
errors in the measurement located a at A and B, s B, such that: 
op 


2. If the standard deviation of the pressure differential, 1 Sane », Obtained from — 


FOA is greater than or equal to the additional pressure differential, A(AP), — 7” 


 —_—- from the length of time the pipe has been in service, then the measurement _ 
devices should have an accuracy capable of providing values in the range of _ 
the mean value of the pressure differential, AP, nk the additional ye 


i 
liferential, O(4/), resulting [rom the period Of time the pipe has Deen in Service, 


1981 


> 


s Consider now how the results obtained in n the preceding from from n the application — 
of FOA can be used to select accuracy needs in pressure measurement devices 
_ for field measurements. If ‘‘perfectly’’ accurate pressure recording gauges were 
_ located at the ends of the pipe then with the mean values of the parameters | 
: being considered and neglecting the model error term, , the analysis of the network | 
will indicate a a value of 17.1 KN/m n * for AP. ‘Furthermore, ‘it can be shown : 


that the of the measured pressure differential ary is given by (1): 


pressures A and B respectively; 
= the mean values of measured Pressures at A and ‘B. ‘When a 
“1% accuracy is considered for each pressure gage, the value of « Cor, obtained 
4 s 2.76 KN/m? and for 2% accuracy, the corresponding value is 5.52 "KN/m?. ol 
_ From the results plotted in Fig. | it follows that pressure gages having 7 
2% accuracy or better, will provide values of pressure differential, AP, within | 


FIG. 1.—Variation of Pressure Differential, AP, for Example Pipe ond 


the uncertainty range calculated by FOA ‘(for the example pipe 
if this accuracy | is not attainable, it is ‘recommended that the ‘field experiment 
haps 
“modif ying the lengths of pipes considered in 1 the field testing program. a 
_ It is noteworthy to mention here that the aforementioned FOA method was : 
aiet to many existing large networks, including the City of Hamilton water 
_ distribution network (Canada). In some cases, the required accuracy for sensible ; 
pint measurements was as high as 0.5% (see Ref. 3 for further par Se th a 


‘The application of FOA in determining the sensitivity ofa particular | sinitesion 


= predetermined uncertainties. In this paper, FOA principles were utilized — 


_ to examine the sensitivity of the pressure differential across a pipe due to 
uncertainties in the mathematical model, diameter, discharge, hydraulic slope 
exponent and Hazen-William’s coefficient. The results obtained were also utilized 


to obtain a required accuracy of measurement pei for providing meaningful = 


| 
| 


wi — TECHNICAL NOTES 


results int field measurement ement of pres: pressure sure dif! f erentials 4 along se sections sof distribution 


77 


ite rib 


In conclusion, this study shows uncertainties in the input data of water 
distribution networks are ‘substantial, and therefore could lead to meaningless 


results in the subsequent analysis re regardless of the ‘efficiency of the co computer | 


- program employed; data collection programs must be very carefully designed. ee. 
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: Discussions may be submitted on ~_ ecsiiniie: paper or technical note published i in any 

Journal or on any paper presented at any Specialty Conference or other meeting, the Proceedings 

of which have been published by ASCE. . Discussion of a paper/technical note is open to. 

anyone who has significant comments or questions regarding the content of the paper / technical 
note. Discussions are accepted for a period of 4 months following the date of publication 
of a paper/technical note and they should be sent to the Manager of Technical and Professional — 

_ Publications, ASCE, 345 East 47th Street, New York, N.Y. 10017. The discussion period may 7 

be extended by a written request fromadiscusser, 

The original and three copies of the Discussion should be submitted on 8-1 /2-in. 220-mm) 
by ll-in. (280-mm) white bond paper, typed double-spaced with wide margins. The length of 

_ a Discussion is restricted to two Journal pages (about four typewritten double-spaced pages 

of manuscript including figures and tables); the editors will delete matter extraneous to the 
subject under discussion. If a Discussion is over two pages long it will be returned for shortening. 

_ All Discussions will be reviewed by the editors and the Division’s or Council’s Publications ‘ 
Committees. In some cases, Discussions will be returned to discussers for rewriting, or they 
‘may be encouraged to submit a paper or technical note rather thana Discussion. = 
_ Standards for Discussions are the same as those for Proceedings Papers. A Discussion is 

subject to rejection if it contains matter readily found elsewhere, advocates special interests, 

is carelessly prepared, controverts established fact, is purely speculative, introduces personalities, 
or is foreign to the purposes of the Society. All Discussions should be written in the third — 
_ person, and the discusser should use the term “the writer’’ when referring to himself. The 
author of the original paper/technical note is referred to as “‘the author.” 
i Discussions have a specific format. The title of the original paper/ ‘technical note appears 
at the top of the first page with a superscript that corresponds to a footnote indicating the 
} | month, year, author(s), and number of the original paper/technical note. The discusser’s — 
| . ‘name should be indicated below the title (see Discussions herein as an example) together with 
ASCE membership grade (if applicable). 

_ The discusser’s title, company affiliation, and business address should appear on the first 
) i. of the manuscript, along with the Proceedings paper number of the original paper / technical 
} “note, the date and name of the Journal in which it appeared, and the original author’s name. a 
P - Note that the discusser’s identification footnote should follow consecutively from the original | 
"Paper /technical note. If the paper / technical | note under discussion contained footnote 


Figures supplied by the uni should be designated by letters, starting with A. This also 
“applies separately to tables and references. In referring to a figure, table, or eens that 
jun in the original paper/technical note use the same number used in the original. a 
a It is suggested that potential discussers request a copy of the ASCE Authors’ Guide to 
the Publications of ASCE for more detailed information on preparation and submission of 
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nportan of four wellknown probabihty 

INDUCED BY TURBULENCE IN STILLING BASIN" 


With the by Rangaswami Narayanan> 


- that the method presented is over-simplified for the obvious reason that there 
S are so many uncertain factors. The uncertainties ‘such as the duration of pressure oe 
‘below vapor pressure, the onset of cavitation coincident at Pp, itself, and the 
others listed in the paper, make a rigrous analysis of the criterion for inception - 
of cavitation extremely difficult. Obviously, further research is necessary before ‘ 
a Predictive method is developed for the inception in stilling basins Goling 7 


_ The writer’s assumption with r respect | to the mean | pressure re distribution in 
a the hydraulic jump is not all that drastic. The discussor’s result regarding the — 
pressure distribution is at variance with the results of Rajaratnam (8). Rajaratnam — 
presents Static pressure distribution within the hydraulic jump to show that 
' near the bed for about 5% of the depth the pressure is equivalent to what 
would be hydrostatic. If at all, the pressure is less than hydrostatic in the 
_ body of the hydraulic jump. ‘The deviation from the hydrostatic distribution 
is concentrated mostly in the earlier part of the jump. But in the region where - 
__ the pressure fluctuations are most active, the assumption of pressure distribution - > 
as assumed in the paper appears sensible, 
- The discussor argues in the main that the intermittency could be smaller 
than that deduced in the paper. While accepting the points raised by the discussor, _ 
the writer “would like to point out that the inception of cavitation as coincident — 
_ with the existence of vapor pressure is conservative. Clear evidence is now — 
available that the inception occurs at pressures much higher than the vapor 
_ pressure at the ambient temperature. Therefore, the intermittency could be larger » 


‘The writer, assuming that the duration of the instantaneous pressures below 
the level p, is typically of the same order as the inverse of the dominant frequency " 
of pressure pulsations, finds that ¢,, could be greater than ¢.,. Thus, it is 


appropriate to assume that there is enough time for the bubbles to gow 


8. Rajaratnam, LN, “The — Jump as a Wall Jet,”’ nimnabe of the » Hydraulics Division, 
ASCE, Vol. 91, No. Proc. Paper 4482, Sept., 1965, pp. 107-132. > 


‘April, 1980, by Rangaswami Narayanan (Proc. Paper 15305). 
*Lect., Dept. of Civ. and Struct. Engrg., Univ. of Manchester inst. of Sci. and Teck. 
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* -—The fol llowing © corrections should be made to the paper: 

ce nee 618, line 4 from bottom of the page: Should read y, = y 
33.77 m; = 5.7 33.77 = 192 m; and y,, = 


Loc PEARSON Tyre 3 DisTRIBUTION: 


Closure by Donthamsetti Veerabhadra Rao,’ M. 
would like to th thank. for his interest in paper and for 
i. his valuable comments. The writer concurs with the discusser’s view that ever 
_ since the Water Resources Council (WRC) recommended (5) the use of LP 
- for flood flow frequency analysis in ee much knowledge has been added = 
in respect of LP andits fitting procedures, = 
i, Until Bobee (1) derived the density function of LP in 1975, the only known 
_ procedure (by moments) for fitting LP \ was through the logarithmic moments _ x 
i _ (LGMO). Let Y, = In X,, in which X, = ith item of real data. Then the 
logarithmic SP, mean (Y), variance (S? ), the coefficient of skewness (CS,) 
were calculated from the log data Y,, and the parameters of LP were estimated « 
a from P. (See LGMO method, Ref. 13.) The efforts of different investigators _ 
were in the direction of improving this LGMO method; it was recognized that es 
a the sample CS, was biased, and weighting C S, on the basis of a generalized f 
skew map was ’ suggested (10). Refs. 1 and 13 "show that alternative methods 


ha of moments are now available for fitting LP to data samples parameters of he 
‘ LP can also be estimated by the moments of real data (RLMO), or by mixing ~ ¥ 

the real and log moments. In the WRC method the coefficient of variation ie 
| in log space is considered accurate and the logarithmic sk: skew is adjusted. Logically, ah 
4 this consideration applies also to data in real space, i.e., the c coefficient of as 


_ variation of real data may be considered accurate. Thus, there are four sample | 
(4 SP which may be regarded as accurate or less biased, viz., X, S?, Y, and tt 
=» Since LP has three parameters to be estimated, it is possible to preserve | 
three of the preceding SP while fitting LP to a ‘sample. This was the concept 
_ a behind the development of the ‘methods of mixed moments (13). The method 
D MXMI referred to by the discusser preserves X, S?, and Y of the tae 


“May, 1980, by Donthamsetti Veerbhadra Rao (Proc. Paper 15391), 
*Hydro., St. Johns Management District, P.O. Box 1429, Fla 


Xa 
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has become (15). In Ref. 15 the writer analyzed and ‘compared some 


_ important properties of four well known three-parameter probability distributions, © 
--Ramely, lognormal, Weibull, P, and LP. In this study LP was found to possess — 
some appealing properties (not found in other distributions) in the context of | '; 


- fitting LP to real data by ML is equivalent to fitting a P distribution to logarithmic 
data by ML. Matalas and Wallis (14) studied at length the ML estimation for 
parameters of P, and their concluding remarks were: f 
—several difficulties related to the properties of the P distribution were ol 
encountered in deriving both moment and ML estimates. Although these - a 
. difficulties may not seriously affect the use of this distribution from the . 
point of view of curve fitting, the difficulties are sufficiently troublesome _ a 


_ that using other distributions warrants consideration, ‘particularly if ML 


Bobee and Robitaille (8) express that for P distribution ML estimates are not » 
_ jointly sufficient in which case the method is not optimal for small samples. — 
‘ TABLE 7.—Mean and Standard D Deviation of 100- -yr Flood Estimates Based on 1,000 


_ The writer has seen some me advantages when calculations are re performed by using 
| dimensionless data (K,) in fitting LP by ML. Parameter, c, has a positive sign 
4 when LP is bounded upwards and negative sign when bounded in the lower = 
- end. When |c| is large the solution is practically equal to that obtained by 
the t two- parameter ‘Such a an easy of LP parameters 


| 
size Method} Mean | ation | Mean/| ation | Mean ation | Mean yoo 
---:20-—s« |: MXM1} 2.306 | 0.419 | 3.079 | 0.626 | 2.458 | 0.643 | 2.374 | 0.640 > 

"he | 2.208 | 0.443 | 4.025 | 2.459 | 2.548 | 0.695 | 2.280 | 0.583 
30s | MXM1| 2.302 | 0.349 | 3.074 | 0.528 | 2.475 | 0.538 | 2.406 | 0.565 
the | 2.217 | 0.364 | 3.298 | 1.459 | 2.599 | 0.622 | 2.358 | 0.519 
: hs sas MXMI | 2.305 | 0.287 | 3.086 | 0.417 | 2.509 | 0.446 | 2.445 | 0.476 | 

2.249 | 0.291 | 2.994 | 0.572 | 2.612 | 0.489 | 2.434 | 0.447 
|. MXM1 | 2.305 0.210 | 3.072 | 0.297 | 2.539 | 0.331 | 2.479 | 0.362 
— ML | 2.282 | 0.213 | 2.995 | 0.274 | 2.606 | 0.339 | 2.486 | 0.328 — 

; 7 described in Ket. 14 the writer obtained M estimates for the P data Sampies 7 
generated by Monte Carlo simulation presented in Ref. 13. For the four parameter — 
| 8 sets used in the preceding simulation experiments Table 7, 


_ Table 7 shows that for parameters set Nos. 1-2 (y, negative) the estimates | 
_ of MXM1 are less biased and less variable than the ML estimates; for small 
- ‘sample sizes (n = 30) the ML estimates | for several of these samples — 
re} be considered unacceptable. When jy, is positive (parameters set Nos. 3-4) 
the MXMI and ML estimates have somewhat balancing attributes. In general, | 
the ML method gives satisfactory estimates only for large samples not commonly 
_ found in practice. Thus, MXM1 method which is simpler in its application and 
gives better estimates compared to the other moments methods as well as ML, 
_ is preferable for LP parameter estimation. Complete details of the writer’ Ss 
- approach to ML method and the Monte Carlo simulation results will be presented — 


14. Matalas, N. C., and Wallis, J. R., “‘Eureka! It Fits a Pearson Type 3 Distribution, a 
wal Water Resources Research, Vol. 9, No. 2, Apr. , 1973, pp. . 281- 289. — a 
15. Rao, D. V., ‘“‘Three-Parameter Probability Distributions,’ Journal of the Hydraulics 2 
Division, ASCE, Vol. 107, No. HY3, Proc. Paper 16124, Mar., 1981, pp. 339-358. 
Page 854, line 1: Should Read “SP” instead of “of SP” enw 
‘Page 859, paragraph 1, line 3: Should Read “‘Note that K »”” in instead of ‘‘Note . 
Page 860, “paragraph: 1, lines 5-6: Should read values which are 
than | ,000-yr K values by LN2 and gamma distributions” = of ‘ “eee 
K,, values, which are less than distributions” 
_ Page 868, paragraph 3, line 7: Should read **K values for a given T; (a) a ie 
(b) increase and then decrease; ‘te instead of “K values for %, (a), — 
(b), increase and then decrease, | | 


Page 871, paragraph 2 (of Summary), line 8: Should read seg or the upper 


quantiles,’ instead of ‘‘For the upper 


CHARACTERISTICS OF 


and Needs of the Committee on Surface-Water Hydrology. 

_ North states that the paper clearly emphasizes the urgency of developing 
alternatives to the M-day, T- -year low flow approach, and he describes a a 


P:! of characterizing low flow. by the duration of negative : daily flow runs. That 


m= 1980, by Task Committee on Low-Flow Evaluation, Methods, and Needs of PY 


the Committee on a Surface- Water ene of the Hydraulics Division (Proc. Paper 15400). — 


| 
BASIL. 


method should be further atihiiaihes so so that i its ts utility can be assessed. _ 
weakness appears to be that few runs are available for defining the rer aere 
curves of runs for low values of threshold flow. pom, 
... A study report by the Institute of Hydrology (43) defines a low flow oe 
_ either as the length of period that the flow is continuously below a threshold 
or as the deficit volume that would be required to maintain the flow at threshold. — 
_ Frequency curves of these are prepared for various threshold discharges. > phot 


a _ Alternatives to the M-day, T-year low flow are needed because no one descriptor 
i of low flow characteristics is suitable for all hydrologic applications. Our paper 
- pointed out, conclusion 4, that hydrologists should determine the most suitable 


ways of using low flow information for various applications and then provide 


a We appreciate North’s contribution to the Task Committee Teport. Shel got 


oni ad Dive 
Reports 1-3, “Institute of Hydrology, Low Flow Studies, “Institute of ‘Hydrology, 


°F. ASCE 


_ Wilson may have misread the problem posed in the paper. The ; analysis w was 


undertaken to deal with cases where conventional boundary layer calculation — é 
methods” are not applicable; some - examples of such flows were cited in the — 
paper. The specifically the discusser’ s 9 and 10 


‘the drop in the specific energy line, then it is necessary to use coefficients — : ~ 


: like those used by Jaeger (10). However, for estimating changes in depth and ste 
. ‘mean velocity in complex developing flows, it is much simpler to guess a value © 
of displacement thickness than it is to guess these other coefficients; this is 
the crux of the argument in the paper. If the discusser believes he can {reat 
complex developing flows by the method he advocates, let him reproduce 
_ Kindsvater’ s data (8) as was done by the writer and shown in the inset of 
1980, by Edward Silberman (Proc. Paper 15516), 


*Prof., St. Anthony Falls Hydraulic Laboratory, University of Minnesota, Minneapolis, _ 


| 

4 

Ly 
the discusser’s claim to the contrary, displacement thickness as given by Eq. : 

. _ 2 is positive definite under the conditions stated on the third line from the _ : 

bottom of page 1237 of the paper. Displacement thickness enters into the definition _ 

_ of critical depth via Eas. | to 4 and not because the writer added it arbitrarily, 

| . 
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‘ Bar or Gravet-Bep STREAMS” 


_ The authors are to be congratulated on tackling the problem of flow resistance 
at low flows in gravel-bed rivers since so often previous work has concentrated — x 
on conditions at bankfull flows. However, while the authors’ approach makes ‘ 

a useful start on the problem the result is, as noted, an empirical formulation ; 

for a bulk bar resistance which is subject to uncertainties, typified by the scatter a 

in Fig. 5 and by the slightly arbitrary criterion that, once sediment movement — ce 


begins, bar resistance is negligible. In order to produce a more generally apy applicable _ 
equation it will eventually be necessary to consider the processes by which 
i bar resistance is generated and this writer would therefore like to comment 
= on the processes most likely tobe of importance. 
Firstly, in addition to the bar and grain resistances there are resistance 
components | related to channel cross-sectional shape, secondary circulation and 
: nonuniform shear stress distribution (10,40). Consequently, it is an approximation 
to assign all the resistance other than grain resistance to the influence of bars, — 
as the authors have done. However, at low flows bar resistance may well be > 
the most important of the ‘nongrain resistances so the may 


- Secondly, it may be conjectured that there are two main processes _ which a 


bars increase the flow resistance. The first is the ponding effect of the bars | 
or riffles which, because the longitudinal bed profile is nonuniform, can act 
as controls on the flow in the pools, keeping those flows deeper and slower 
than would otherwise be the case. The second resistance is the result of the 
shallow flow over the bars or riffles which is characterized by features of 
large- -scale | roughness : and an increased resistance to flow. As s pools are usually . 
longer than riffles and as ponding usually retards the flow more than does 
large-scale roughness, the contribution of the riffle flow to the overall resistance 

of a river reach is less than that of the ponding effect in the pools, although 
of course the grain resistance of the riffles is higher than that of the pools, — 
Both resistance effects diminish as depth | increases. aA 
_ Means of accounting for large- -scale roughness are discussed elsewhere (39) 

. so only the ponding effect is considered here. It is proposed that this be quantified | 7 
using the residual depth, or the depth in a pool when discharge is zero. Because = 
this residual depth exists, the depth at a given discharge is greater than it would — 
be for the same discharge over the same bed material with a uniform bed ; 
"profile, while the mean velocity is lower. ‘However, as ‘discharge | and depth 
increase, “the ra ratio of residual depth to actual depth decreases and consequently 
the relative effect of the residual depth (i.e., the degree of ponding) decreases. eS 
lc possibility of using residual depth to account for the major portion of | 
October, 1980, by Gary Parker and Allan W. Peterson (Proc. Paper 15733), 

"Sr. Scientific Officer, Inst. of Hydrology, Wallingford, Oxon, United Kingdom. oy 4 


| 
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_ bar resistance is vis Miaieteaied here u using data from om the uy upper River Swale: in England 

G8). The data were collected from a reach of channel 1,100 m (3,609 ft) long © 

a 8 pool/riffle units) at three discharges. Average value of D,, for 3 
é the bed material [based on 19 Wolman-style samples (41)] is 132 mm ©. 433 ae 


ft), : average channel slope is 0.0027 and channel width is in the range 17 m-20 m 
(55.8 ft-65.6 ft). The reach-averaged residual hydraulic radius J, (or hydraulic 
radius at zero discharge) was obtained by back extrapolation on a linear plot a 
of average reach hydraulic radius J on discharge Q, giving the value of - 
d (=0.2 m or 0.656 ft) at which Q = zero. This method is rather approximate he 
given only three data points but the value so derived is suitable for 
; Following the authors assumption that flow resistance in the absence of bars 
6 consists mainly of grain resistance, hypothetical values of discharge Q,, were 
i. calculated for the three flows as if only grain resistance were present. For 
this calculation the Hey equation (10), rather than Eq. 2, was used since it 


has given good results in the fiel field: LW 


in i A= = flow cross-sectional area; g = ‘seaiaiii due to gravity; S 
op cecigin vines of D ip tee 
E 1.—Average Flow Parameters for River Swale Field 


discharge, Q,| discharge, -Reach-averaged 
incubic | in cubic hydraulic 

meters per meters per radius, |, 


| 


= 
Note: 1m = 3.281 ft; 1m’ 


= channel ‘slope; and a = a channel shape parameter. If bar resistance were 


_ negligible it would be expected that Q and Q,, would be approximately equal, 
so the ratio Q@/Q,, was calculated as a rough (inverse) measure of the effect 
of bar resistance. Also, if bar resistance does depend largely on residual —e Mg 
it would be expected that the magnitude of the resistance would be an inverse — 

_ function of the ratio (J — J,)/J. Then, assuming that other resistances are 

> negligible the two ratios (given in Table 1) should vary in a related manner, 

3 fz not necessarily at the same rate. Thus, with Q = 0, and therefore 

» both ratios are zero. _As discharge fi rises, the bars become drowned 


1 
a q a 1 support this hypothesis as far as they go but, as the data are few, the ow 
influence of other reistances cannot be gauged. Nevertheless there seems to a 


e could be calculated directly from the a 
residual depth (the source resistance), which would be more 
satisfactory than deriving it from shear stress criteria which are less —— ; 


q related to the process involved. Considerably more experimental data, though, — - 

will be needed to bring this idea to fruition, alin 
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Discussion Mario Letelier S.,” and Hans J. Leutheusser,* ASCE 


The authors address a problem which has already ¢ 


in the recent technical literature (e.g., 7). In particular, they propose an 
_ approximate method of analysis for ey laminar closed- conduit flow based 


‘the flow is at any instant subject to the corresponding steady- -state friction. 
_ In order to properly assess the validity of the authors’ equation of motion, | 


“viz. Eq. 5, it is convenient to write it in nondimensional form, i.e. | 


which Q* = Q/Qo; = Ap* = Ap/pLv and and 
T, are reference values of Q and ¢, 


- os 1980, by David Stavitsky and Enzo Macagno (Proc. Paper 15922), 

Prof., Departamento de Mecanica, Facultad de Ingenieria, Universidad Técnica del 

 “Prof., Dept. of Mech. Engrg., Univ. of Toronto, Toronto, Ontario, io, MSS 1A4, Canada. 
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DISCUSSION 
is an “‘unsteadiness aumber,” which can, , for instance be unin’ as th 


between some “‘pipe-t “fluid” time 


, iLe., R?/v, and an external ‘‘forcing’’ 


tom time :.. Clearly, when 2 — 0, Eq. 1 reduces to its steady-state equivalent. iN 
_ A more general and exact form of the one-dimensional momentum equation 
¥ than Eq. 32 can be deduced from the Navier-Stokes equation. Thus, for yoinu sil 


flow of an incompressible Newtonian fluid, there ae ates ene 


+ — 


ose 
which = u(r,t) is the local axial velocity which is a function of the radial 


—-» 


_ After some manipulations (S, 6), Eq. 34 can be transformed = as series 
expression, viz. 3% 


+ 
which i is valid for any - value of { a. A somewhat simpler form of 35 which, — 
however, converges only to certain values of ai is the following: 


The authors’ model, Eq. 32 of this is a version 
of 36 with an incorrect coefficient for the first time-derivative of Q*. In the = 
Pe light of Eqs. 35 or 36, this model can be « expected to yield acceptable results 
only for small values of 0, i.e., when all terms of second and higher order 
i in series Eqs. 35 and 36 are negligibly small. One interesting exception > 
for = o. In this case, transforms into the inviscid (i. e., ideal flow) 


form of Eq. 34, viz. 


dQ 


__When applied t to the case e “Starting Flow i in Pipe,” the authors mode, Ea 


using only the first two terms on the left-hand side, then any error would 

- be negligible. The reason for this behavior derives from the fact that a “starting 

flow”’ involves large acceleration at small time and, therefore, many time- -deriva- a 

tives are necessary to correctly describe the motion at its beginning. 

_ In the second example given in the paper, viz. ‘Oscillatory Flow in Pipe,” 

‘the authors’ model reproduces exact results for both 0 and », i. — 
“‘steady state’’ and “‘ideal flow,” respectively. Nevertheless, for other values 
- of 2, Eq. 32 does entail considecabie error. Thus, according to computations _ 


by the Eq. 32 yields, e.g., for = = an which 


9 

aia 

e 
| 
d 
7 
4 
rate of flow, which is in error by relative to the exact time (6). On the 


rise (8), show that the artnet of the authors’ model will invariable entail 
_ nonnegligible errors. Nevertheless, the discussers readily concede that equations — 
of the form of Eq. 32 are very convenient for obtaining a first, rough assessment 


of a particular transient flow problem. — ‘However, and _ especially when large © 
accelerations exist in the flow, more precise computations call for the use of © 
"equations such as Eq. 35.0r36, 
> 
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- is in error by 17% relative to the exact value. And this error can be still larger £ 
7 for other values of 9 (5). On the other hand, Eq. 35 reproduces the exact 2 
| e solution for all values of 2, while Eq. 36 reproduces it up to 2 = ee 
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